% Résumé du stage de M2 de Rida sur le défilement des tiges de Douglas

\documentclass[12pt,french]{report}

%

\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}

\usepackage{url}

\usepackage{graphicx}

\usepackage{natbib}

\usepackage[francais]{babel}
%\usepackage[french]{babel}

\usepackage{amsmath}
\usepackage{amssymb}

\usepackage{rotating}

\usepackage{tocbibind}

\usepackage[firstpage]{draftwatermark}

\usepackage{pdfpages}

\usepackage{placeins}	% pour \FloatBarrier
%\usepackage{empheq}

\graphicspath{{Fig/}, {fig/}}

\selectlanguage{french}

%document principal

\SetWatermarkText{Version provisoire en cours de modification}
\SetWatermarkLightness{0.5}

\begin{document}

\title{Qualité du Douglas en France}

\author{Fleur Longuetaud et Frédéric Mothe 
\\avec la participation de Daniel Rittié}

\date{Septembre 2014}

\maketitle

\newpage

\tableofcontents

\newpage

\chapter{Présentation des jeux de données}

\section{Jeu 1 : Ecouves}

\label{Ecouves}

\subsection{Présentation générale}

Ce jeu de données est constitué de 27 Douglas de la forêt d'Ecouves (61) (coordonnées Lambert II étendu : 429 400 E  ; 2 425 040 N). Il s'agit d'une plantation (2 m $\times$ 2 m ; 2500 tiges/ha ; plants de 3 ans plantés pendant l'hiver 1950/1951) et tous les arbres ont 63 ans au moment des mesures (âge fin 2010 ; arbres coupés en 2011). Le dispositif a été installé alors que les arbres avaient 18 ans en 1965.

Trois traitements sont représentés : (A) un témoin (9 arbres numérotés de 10 à 18), (C) un traitement correspondant à des éclaircies fortes au profit des arbres d'avenir (9 arbres numérotés de 19 à 27) et (B) un traitement intermédiaire de type ONF local (9 arbres numérotés de 1 à 9). 

Pour chacun des traitements trois statut sociaux sont représentés soient 3 arbres dominés, 3 co-dominants et 3 dominants.

Pour chacun des 27 arbres, une analyse de tige détaillée a été réalisée avec mesures des largeurs de cernes et de la limite aubier/duramen sur quatre rayons à plusieurs hauteurs.

Les autres variables disponibles sont : le $C_{130}$, la hauteur totale, la hauteur de la base du houppier, la hauteur de la 1ère branche verte, la hauteur du 1$^{er}$ frottis, la projection du houppier, les positions et diamètres des voisins.

Les éclaircies ont été réalisées en 1969, 1975, 1980, 1985, 1992 et 2000.

Note 1 : pour les mesures de surface terrière, de $H_0$ et les effectifs, les données pour un traitement sont obtenues à partir de la moyenne de trois placettes du traitement.

Note 2 : le peuplement témoin non éclairci sera considéré comme hyperdense et utilisé pour le calcul de la droite d'auto-éclaircie même si des chablis ont eu lieu notamment en 1999.

Note 3 : à première vue les arbres ont souffert en 1976 mais pas en 2003.

\subsection{Caractérisation des peuplements de l'étude}

\subsubsection{Droite d'auto-éclaircie}

On estime dans un premier temps la droite d'auto-éclaircie (Équation (\ref{Reineke}) ; Fig. \ref{auto_eclaircie}) à partir des trois placettes du peuplement A non éclairci (Reineke, 1933 in Pretzsch, 2009). La droite obtenue par Reineke (1933) sur des peuplements de Douglas équiennes non éclaircis aux USA (états de Washington et Oregon) a été tracée (en noir). Un ajustement aux moindres carrés (OLS) a été testé (en supprimant les premières années car le peuplement planté n'était pas hyperdense ; en vert). Le meilleur ajustement (validation visuelle ; en rouge) est donné par la méthode de l'analyse de frontières stochastiques (SFA ; package \textit{frontier} sous R).

\begin{equation} 
	\log(N) = 10.90 - 1.15 \log(D_{g})
	\label{Reineke}
\end{equation}
avec $N$ le nombre de tiges (tiges/ha) et $D_{g}$ le diamètre quadratique moyen (cm). Les valeurs pour l'intercept et la pente sont celles obtenues par la méthode SFA.

\begin{figure}[!h]
	\centering	
	\includegraphics[width=9cm]{auto_eclaircie.png}
	\caption{Ajustement de la droite d'auto-éclaircie à partir du peuplement A.}
	\label{auto_eclaircie}
\end{figure}

\subsubsection{Trajectoires de $RDI$}
La droite d'auto-éclaircie ajustée précédemment nous sert ensuite à tracer les trajectoires de $RDI$ (Equation (\ref{equa_RDI})) des trois peuplements (Fig. \ref{RDI}). La différence de $RDI$ entre les deux traitements éclaircis est de l'ordre de 0.1. L'effet de la tempête de 1999 est visible sur les trajectoires de RDI. 

\begin{equation} 
	RDI = N / N_{max} = N \times D_{g}^{1.15} / \exp(10.90) 
	\label{equa_RDI}
\end{equation}

\begin{figure}[!h]
	\centering	
	\includegraphics[width=9cm]{trajectoires_RDI.png}
	\caption{Trajectoires de RDI pour les trois traitements A, B et C (trois placettes par traitement).}
	\label{RDI}
\end{figure}

\subsubsection{Croissance en hauteur dominante}
La croissance en hauteur dominante des peuplements est représentée Fig. \ref{H0}. On peut considérer que les trois peuplements sont de fertilité comparable. En superposant des courbes de croissance en hauteur dominante ($H_{0}$) de référence pour le Douglas \citep{Angelier2007}, nous pouvons estimer que la fertilité des peuplements est de classe 2 ($H_{0}$ à 50 ans = 33 m).

La première éclaircie a eu lieu en-dessous de 15 m de $H_{0}$ pour les deux peuplements éclaircis alors que pour une fertilité de classe 2 une première éclaircie vers 17-18 m est préconisée avec des rotations de 7 ans tout au long de la vie du peuplement \citep{Sardin2012}.

\begin{figure}[!h]
	\centering	
	\includegraphics[width=9cm]{croissance_H0.png}
	\caption{Croissances en hauteur dominante ($H_{0}$) pour les trois traitements A, B et C (trois placettes par traitement). Les courbes de croissance en hauteur dominante de référence (classes 1, 2 et 3 du guide ONF) sont superposées. Les lignes grises verticales indiquent les âges pour lesquels des éclaircies ont eu lieu.}
	\label{H0}
\end{figure}

\subsubsection{Comparaison au référentiel des sylvicultures ONF}
Page 16 du correctif du guide de sylviculture 2012 \citep{Sardin2012}, deux graphiques sont donnés représentant respectivement l'évolution de la surface terrière ($G$) et du nombre de tiges par hectare ($N$) après coupe en fonction de $H_{0}$ pour la classe de fertilité 2. La Fig. \ref{G_N_vs_H0} a été obtenue à partir de nos trois peuplements. Les trajectoires sont très différentes entre le référentiel ONF et les mesures de terrain réalisées sur les deux peuplements éclaircis. Pour $H_{0}$ = 40 m, le référentiel donne un $G$ de l'ordre de 31 m$^{2}$/ha alors que les peuplements B et C pour la même valeur de $H_{0}$ sont plutôt légèrement au-dessus de 60 m$^{2}$/ha. Le nombre de tiges est également bien plus élevé que ce qui est recommandé. Nos peuplements ne sont donc a priori pas très représentatifs des sylvicultures préconisées dans le guide ONF. Ces constatations sur la forte productivité du dispositif d'Ecouves ont également été faites par d'autres collègues du laboratoire.

\begin{figure}[!h]
	\centering	
	\includegraphics[width=6.5cm]{G_vs_H0.png}
	\includegraphics[width=6.5cm]{N_vs_H0.png}
	\caption{Surface terrière après coupe (m$^{2}$/ha) en fonction de $H_{0}$ (m) (à gauche). Nombre de tiges par hectare après coupe en fonction de $H_{0}$ (m) (à droite).}
	\label{G_N_vs_H0}
\end{figure}

\FloatBarrier

\section{Jeu 2 : Carré latin d'Amance}

\subsection{Présentation générale}

Ce jeu de données est constitué de 24 Douglas provenant d'une plantation expérimentale de la forêt domaniale d'Amance (54) communément appelée "Carré Latin d'Amance" \citep{Oswald1976} (coordonnées Lambert II étendu : 895 970 E ; 2 425 040 N). La plantation comprend 16 carrés de 10 ares avec différentes modalités de densités initiales et d'éclaircies.  Tous les arbres ont 39 ans au moment des mesures (âge fin 1991 ; arbres coupés en 1992). Le dispositif a été installé alors que les arbres avaient 3 ans en 1955.

Quatre traitements sont représentés par 6 arbres chacun : 
\begin{itemize}
\item (B1) planté à 2 m $\times$ 2 m ; 2500 tiges/ha, non éclairci,
\item (C3) planté à 2,50 m $\times$ 2,50 m ; 1600 tiges/ha, une éclaircie en 1973,
\item (D1) planté à 3 m $\times$ 3 m ; 1111 tiges/ha, non éclairci,
\item (D4) planté à 3 m $\times$ 3 m ; 1111 tiges/ha, trois éclaircies en 1973, 1979, 1983.
\end{itemize}

Les 6 arbres de chaque modalité ont été choisis en considérant leur hauteur de base du houppier (HBH) : 2 arbres de HBH proche de la moyenne du peuplement, 2 arbres de HBH proche de la moyenne plus un écart-type et 2 arbres de HBH proche de la moyenne moins un écart-type.

Pour chacun des 24 arbres, une analyse de tige détaillée a été réalisée avec mesures des largeurs de cernes sur quatre rayons à plusieurs hauteurs.

Les autres variables disponibles sont : le $C_{130}$, la hauteur totale, la hauteur de la base du houppier (\textit{"lowest whorl wearing about 75\% of its nominal needle mass"}), les positions et diamètres des voisins ainsi que les projections de houppiers (sauf pour la modalité D4) \footnote{Les données concernant les voisins et les projections de houppiers n'ont pas encore été retrouvées.}.
De plus, des mesures d'infradensité ont été effectuées par groupes de 5 cernes sur quatre rayons de la plupart des rondelles prélevées.

\chapter{Comprendre et modéliser le défilement des tiges de Douglas en France}

Ce travail a été initié lors du stage de M2 de Rida Kheirallah financé par la convention Modelfor 2012-2015 puis approfondi par la suite.

\section{Revue de la littérature}

\subsection{Les différents modèles de profils de tiges}

Pour décrire la forme d'une tige on peut utiliser un coefficient de forme (F) avec $V = FGH$ où $V$ est le volume, $G$ est la surface terrière et $H$ est la hauteur de l'arbre \citep{Fonweban1997} ou bien une fonction de défilement appelée profil de tige.

Un profil de tige est une équation du type : $d = f(h)$ ou $s = f(h) $, où $d$ et $s$ sont respectivement le diamètre et la surface à une hauteur $h$. 

Les fonctions dépendent généralement du diamètre à 1,30 m ($D$) et de la hauteur totale de l'arbre ($H$) : $d = f(D,H,h)$ ou $s = f(D,H,h)$. 

Note 1 : Dans certains cas on peut avoir besoin d'inverser un profil de tige afin d'obtenir une relation du type $h = g(d)$ (cf. ANR EMERGE).

Note 2 : On pourra également lire le document de l'IFN sur les modèles de profils de tiges \citep{Delord1984}.\\ 

Il existe dans la littérature scientifique un grand nombre d'équations de profils de tiges. 

Dans \cite{Fonweban1997} on trouve une revue assez complète des formes d'équations utilisées : de simples fonction de défilement \citep{Kozak1969,Omerod1973} ; des polynômes \citep{Bruce1968,Goulding1976} ; des polynômes segmentés pour tenir compte de la variation de forme au niveau de la souche \citep{Max1976,Demaerschalk1977}. Un certains nombre d'inconvénients de ces fonctions ont été relevés, notamment l'intégration de ces fonctions pour le calcul de volumes n'est pas toujours aisée. 

Toujours d'après \cite{Fonweban1997}, les fonctions à forme variable (\textit{variable-form taper functions}) \citep{Kozak1988,Newnham1988} sont les innovations les plus récentes et l'approche qu'ils ont retenu dans leur étude. Ce sont des fonctions du type : 
\begin{equation}
d(h) = D \cdot x(h)^{\theta(h)}
\label{equa2}
\end{equation}
où $x(h) = (H-h)/(H-1,3)$ et $\theta(h)$ est un exposant qui varie continûment le long de la tige, de la souche jusqu'à l'apex, et qui rend compte de la forme locale de la tige : néloïdique ($\theta(h)>1$), conique ($\theta(h)=1$) ou paraboloïde  ($\theta(h)<1$) \citep{Fonweban1997,Valentine2001} \ref{schema}.

\begin{figure}[!h]
	\centering	
	\includegraphics[width=9cm]{schema.png}
	\caption{Représentation schématique d'une fonction de profile à forme variable.}
	\label{schema}
\end{figure}

Ces fonctions permettent de bien modéliser les variations de forme à l'empattement avec peu de paramètres par rapport aux fonctions polynomiales. De plus avec ces fonctions le problème de la recherche des points d'inflexion ne se pose pas, contrairement aux modèles segmentés. Ces fonctions peuvent néanmoins poser des problèmes d'intégration pour le calcul des volumes : le terme aléatoire d'erreur est négligé et l'ajustement est réalisé sur le diamètre et pas sur le volume. Une solution serait d'ajuster une équation de défilement de la surface terrière et à employer des modèles d'erreur tenant compte de l'auto-corrélation le long de la tige \citep{Fonweban1997}. Pour les équations proposées dans \cite{Fonweban1997}, l'intégration explicite n'est pas possible. 

Pour l'ajustement de l'équation (\ref{equa2}), ils obtiennent $\theta(h) = \exp(0,3409 + 0,1210 / h + 0,1805 (D/H)^2)$. La variable $D/H$ permet de prendre en compte le statut social de l'arbre et la variable $1/h$ permet d'améliorer l'ajustement du modèle au niveau de l'empattement \citep{Newnham1991}. 

\cite{Li2010} listent 10 équations plus ou moins complexes issues de la littérature dont \cite{Max1976}, \cite{Kozak2004}, \cite{Bi2000}, \cite{Valentine2001} et \cite{Sharma2004}.

Dans ce paragraphe, nous décrivons le modèle de \cite{Valentine2001} qui semble intéressant car c'est un raffinement des fonctions à forme variable. Trois segments géométriques sont considérés : un segment basal, un segment central et un segment apical. L'utilisation de commutateurs numériques (\textit{numerical switches}) permettant de modifier la valeur des exposants permet une transition graduelle ou abrupte entre les segments. L'hypothèse est faite que la transition entre les deux derniers segments correspond à la hauteur de la base du houppier. La corrélation des résidus est prise en compte. Le $1^{er}$ segment a une forme neloïdique. Le segment central a une forme de paraboloïde. Le dernier segment a une forme de cône. Une revue rapide de la littérature est présentée dans cet article, assez similaire à celle de \cite{Fonweban1997}. Ils proposent un modèle à forme et à exposant variables dans lequel les paramètres sont faciles à interpréter. C'est la surface transversale qui est modélisée (mieux pour le calcul des erreurs quand on passe au volume par intégration) :
\begin{equation}
A(h) = A(1,37) \cdot (\frac{H-h}{H-1,37})^{\alpha_{1}+S_{1}(h)} \cdot (\frac{H-h}{H-HCB})^{\alpha_{2}S_{2}(h)}
\end{equation}  
où
\begin{equation}
S_{1}(h) = \frac{\theta_{1}H}{1+[h/(\theta_{3}H)]^{\lambda_{1}}}           
\end{equation}
pour $0 \leq \theta_{3} \leq HCB/H$,
et
\begin{equation}
S_{2}(h) = \frac{(h/HCB)^{\lambda_{2}}}{1+(h/HCB)^{\lambda_{2}}}
\end{equation}

\subsection{Les effets des traitements sylvicoles sur la forme des tiges}

Là encore la bibliographie sur le sujet est assez vaste. Nous avons regroupé dans un tableau plus d'une trentaine de références concernant essentiellement des essences résineuses et la liste n'est pas exhaustive. 

Certaines références ne s'intéressent qu'au défilement de la tige et dans ce cas le rapport $H/D$ ou $D/H$ est largement utilisé.

D'autres travaux vont plus loin en ajustant des équations de profils de tiges (le plus souvent des équations de forme à exposant variable, \textit{variable-form taper functions}) et en testant statistiquement l'effet de variables reflétant le traitement sylvicole sur certains paramètres du modèle.

Les variables sylvicoles ayant fait l'objet d'études concernant leur impact sur le défilement et la forme des tiges sont : la densité de plantation, les éclaircies pré-commerciales \footnote{pratiquées à l'est en Amérique du Nord pour les peuplements issus de régénération naturelle}, les éclaircies réalisées avant la fermeture du couvert (\textit{respacing}), l'intensité des éclaircies, la fertilisation à l'azote, le mélange (\textit{nursing mixture}) et l'élagage.
	
Toutes les études confirment une augmentation du défilement des tiges et donc une qualité de forme dégradée quand l'espacement entre les tiges et/ou l'intensité des éclaircies augmentent. Un changement d'allocation le long de la tige est observé quand la croissance est favorisée avec une allocation plus forte dans le bas des tiges que dans le haut \citep{Valinger1992,Mitchell2000,Goudiaby2012}. Certaines études ont montré que la fertilisation n'avait pas d'effet sur la forme des tiges \citep{Thomson1984,Weiskittel2006}. L'effet de la fertilisation, quand il est significatif, tendrait plutôt à augmenter le défilement \citep{Karlsson2000}. Par ailleurs, \cite{Valinger1992} a montré sur \textit{Pinus sylvestris} que la fertilisation favorisait la croissance dans la partie supérieure du tronc. Il y a peu d'études sur l'effet de l'élagage et les résultats ne sont pas très marqués non plus avec soit une absence d'effet \citep{Weiskittel2006}, soit une légère réduction du défilement pour certaines essences (\citep{Henman1963} in \citep{Macdonald2002}). 

Les études concernant le Douglas sont relativement peu nombreuses. Comme pour beaucoup d'autres essences, l'effet des éclaircies est bien marqué \citep{Brix1983,Brix1992,Weiskittel2006} avec le passage d'une forme paraboloïde à une forme plus conique \citep{Weiskittel2006} et une allocation plus forte à la base de la tige, d'autant plus marquée que les arbres étaient petits et élancés au départ \citep{Mitchell2000,Thomson1984}. Un effet de la fertilisation sur la forme des tiges n'a pu être mis en évidence \citep{Brix1983,Thomson1984,Weiskittel2006}. L'effet de l'élagage est soit inexistant \citep{Weiskittel2006}, soit léger avec une tendance à une réduction du défilement (\citep{Henman1963} in \citep{Macdonald2002}).

Les travaux utilisant des analyses de tiges détaillées (\textit{i.e.}, largeurs de cernes à plusieurs hauteurs) restent rares, d'autant plus rares que le nombre de cernes mesurés est important. Par ailleurs, la plupart des études sont relatives à la forme des arbres à un instant donné après une éclaircie ou un régime d'éclaircies mais très peu s'intéressent à l'acquisition de cette forme au cours du temps et à mettre en avant les périodes les plus sensibles concernant le développement de la forme finale. 


\section{Résultats}

Les résultats présentés ci-dessous sont pour l'instant uniquement basés sur les données du jeu 1.

\subsection{$H/D$}

\textit{Dans la version ultérieure de ce document les auto-corrélations seront prises en compte. Cela dit il y a assez peu de chance pour que ça change fondamentalement les résultats.}\\

Le $H/D$ est une variable assez synthétique qui reflète la forme de la tige (conique vs. élancée) et notamment sa capacité de résistance aux vents. De façon un peu similaire à ce qui avait été fait dans \cite{Mitchell2000} nous avons regardé l'impact du statut social et des éclaircies sur les trajectoires de $H/D$. 

Les arbres co-dominants et dominés n'ont pas des trajectoires de $H/D$ significativement différentes l'une de l'autre (Fig. \ref{H_sur_D}). Pour ces arbres, les trajectoires montrent une forte augmentation dès l'âge de 15 ans environ puis une stabilisation vers l'âge de 35 ans environ. En revanche les arbres dominants ont un $H/D$ significativement différent avec une valeur qui reste plus ou moins constante entre 80 et 90 et qui est significativement plus faible que pour les deux autres statuts sociaux indiquant des arbres plus coniques et moins élancés.

Concernant les différents traitements A, B et C, les trajectoires de $H/D$ ne sont pas significativement différentes l'une de l'autre même si le traitement C à fortes éclaircies tend à être un peu en-dessous ce qui semble logique puisque ces arbres sont censés être plus coniques avec une croissance radiale favorisée par les éclaircies (Fig. \ref{H_sur_D}).    

Dans tous les cas il est intéressant de noter que le $H/D$ semble se stabiliser à partir des années 80, c'est-à-dire autour de 35 ans. Est-ce qu'on peut en conclure que la forme de l'arbre s'acquière dans le jeune âge et n'est plus influencée ensuite ?

Pour compléter les observations graphiques, nous avons réalisé une analyse de covariance (ANCOVA) avec ajustement sur l'âge, la modalité d'éclaircies et le statut social ainsi que les interactions entre ces variables. Au final tous les facteurs ressortent significatifs ainsi que certaines interactions. Le statut social arrive en tête (p-value $= 2.2\times10^{-16}$) suivi de la modalité d'éclaircies (p-value $= 8.404\times10^{-5}$). Seule la modalité "éclaircies fortes" (traitement C) est significativement différente du témoin (traitement A) en ce qui concerne la valeur du $H/D$. Un test de Tukey (comparaisons multiples) confirme que C (éclaircies fortes) est significativement différent de A (témoin) et de B (éclaircies modérées) mais qu'il n'est pas possible de trouver une différence significative entre A et B \ref{Tukey_traitement_H_sur_D}. Les différences de $H/D$ entre les traitements C et A et entre les traitements C et B sont de l'ordre de -9.  

\begin{figure}[!h]
	\centering	
	\includegraphics[width=6.5cm]{trajectoire_H_sur_D_par_statut_zoom_bis.pdf}
	\includegraphics[width=6.5cm]{trajectoire_H_sur_D_par_traitement_zoom_bis.pdf}
	\caption{Trajectoires de $H/D$ pour les trois statuts sociaux avec les intervalles de confiance à 95\% autour des moyennes (9 arbres par statut social et par année) (à gauche). Trajectoires de $H/D$ pour les trois traitements A, B et C (trois placettes par traitement) avec les intervalles de confiance à 95\% autour des moyennes (9 arbres par traitement et par année) (à droite).}
	\label{H_sur_D}
\end{figure}


\begin{figure}[!h]
	\centering	
	\includegraphics[width=9cm]{effet_traitement_sur_HsurD.png}
	\caption{Comparaisons multiples pour l'effet de la modalité d'éclaircies sur la valeur moyenne de $H/D$.}
	\label{Tukey_traitement_H_sur_D}
\end{figure}

\subsection{Rendement matière}

\textit{Dans la version ultérieure de ce document les auto-corrélations seront prises en compte. Cela dit il y a assez peu de chance pour que ça change fondamentalement les résultats.}\\

Nous avons calculé le rendement matière comme étant le rapport entre le volume des sciages (\textit{i.e.}, parallélépipède de plus gros volume inscrit dans le billon) et le volume du billon \citep{ONF2013}. Nous avons considéré une bille de pied de 6 m de longueur ce qui semble être classique pour les industriels du Douglas.

Concernant le statut social, les rendements sont les plus élevés pour les co-dominants et les moins élevés pour les dominés.

Les résultats vont dans le sens attendu avec des rendements à terme meilleurs pour le peuplement témoin (A) que pour les peuplements éclaircies (B et C). Cependant, le peuplement B a atteint un fort rendement dans les années 80, meilleur que pour le peuplement témoin, avant de redescendre. 

Comme pour le $H/D$ il y a une stabilisation autour de 35 ans (\textit{i.e.}, dans les années 80). L'âge des arbres semble avoir une influence sur le rendement jusque 35 ans environ et ensuite les différences de rendement observées sont probablement plutôt liées au statut social et au traitement sylvicole.  

Une ANCOVA a montré que seuls l'âge  (p-value $< 2 \times10^{-16}$) et le traitement sylvicole (p-value $= 0.0375$) avaient un effet significatif sur le rendement. En moyenne le rendement du traitement C est de 1\% inférieur à ceux des traitements A et B (le test de comparaisons multiples de Tukey n'est cependant pas significatif). Il convient donc de relativiser car les différences observées à terme ne sont que de l'ordre quelques pour cent.   


\begin{figure}[!h]
	\centering	
	\includegraphics[width=12cm]{Schema_rendement_matiere.png}
	\caption{Schéma d'explication du rendement matière. Source : \cite{ONF2013}.}
	\label{Schema}
\end{figure}

\begin{figure}[!h]
	\centering	
	\includegraphics[width=6.5cm]{Traj_rendement_par_statut_zoom}
	\includegraphics[width=6.5cm]{Traj_rendement_par_traitement_zoom}
	\caption{Trajectoires de rendement matière pour les trois statuts sociaux avec les intervalles de confiance à 95\% autour des moyennes (9 arbres par statut social et par année) (à gauche). Trajectoires de rendement matière pour les trois traitements A, B et C (trois placettes par traitement) avec les intervalles de confiance à 95\% autour des moyennes (9 arbres par traitement et par année) (à droite).}
	\label{Rendement}
\end{figure}

\subsection{Exposant de forme variable $\theta(h)$}

\textit{Dans la version ultérieure de ce document les auto-corrélations seront prises en compte. Cela dit il y a assez peu de chance pour que ça change fondamentalement les résultats.}\\

Dans cette partie on s'intéresse à l'exposant de forme qui permet d'analyser plus finement la forme du profil de la tige indépendamment du rapport $H/D$ puisqu'il met en relation un diamètre relatif $d(h)/D$ avec une hauteur relative $x(h)$.

Partant de l'équation (\ref{equa2}), il est possible de calculer une valeur de $\theta$ pour chaque profil de tige annuel dans l'arbre et pour chaque hauteur $h$ (sauf pour $h$ = 1,30 m et pour $h$ = $H$ où $H$ est la hauteur totale de l'arbre) (Equation (\ref{theta})). \cite{Tasissa1998} utilise la valeur de cet exposant de forme pour étudier l'effet des éclaircies sur la forme des tiges. 

\begin{equation}
\theta(h)=\frac{\log(d(h)/D)}{\log(x(h))}
\label{theta}
\end{equation}

Une ANCOVA avec ajustement sur la hauteur relative $x(h)$, l'âge, la modalité d'éclaircies, le statut social et les interactions deux à deux entre les trois dernières variables a montré que tous les facteurs ressortaient très significatifs sauf l'interaction entre l'âge et le statut social. En particulier, l'effet du traitement est significatif (p-value $= 2.815\times10^{-6}$), ainsi que les interactions entre le traitement et le statut social (p-value $< 2.2\times10^{-16}$) et entre le traitement et l'âge (p-value $= 5.239\times10^{-12}$).

Un test de Tukey (comparaisons multiples) confirme effectivement que les traitements B et C sont significativement différents du traitement A (avec une différence plus marquée entre C et A) mais qu'aucune différence significative n'a pu être mise en évidence entre les traitements B et C (Fig. \ref{Tukey_traitement_theta}).

\begin{figure}[!h]
	\centering	
	\includegraphics[width=9cm]{effet_traitement_sur_theta.png}
	\caption{Comparaisons multiples pour l'effet de la modalité d'éclaircies sur la valeur moyenne de $\theta$.}
	\label{Tukey_traitement_theta}
\end{figure}

La Fig. \ref{Tukey_traitement_theta_HR} illustre les p-values obtenues pour les comparaisons de $\theta$ par classe de hauteurs relatives $x(h)$. Le traitement B (éclaircies modérées) n'est pas toujours significativement différent du traitement A (témoin). En revanche le traitement C (éclaircies fortes) est toujours significativement différent du traitement A, quelle que soit la hauteur relative considérée. Les traitements B et C sont en général significativement différent sauf à l'empattement, c'est-à-dire sous 1,30 m.  

\begin{figure}[!h]
	\centering	
	\includegraphics[width=9cm]{effet_traitement_sur_theta_vs_hauteur_relative.png}
	\caption{P-values pour les comparaisons multiples de l'effet de la modalité d'éclaircies sur la valeur moyenne de $\theta$ par classe de hauteurs relatives. Le seuil de significativité est indiqué pour une ligne en pointillés grise.}
	\label{Tukey_traitement_theta_HR}
\end{figure}

\begin{table}[h]
	\caption{Valeurs de $\theta$ par classe de hauteurs relatives $x(h)$.}
	\centering 
		\begin{minipage}[c]{\textwidth}
		\begin{tabular}{rrrr} 
			\hline
			Hauteur relative\footnote{Les valeurs proches de 0 correspondent à la partie haute de l'arbre et les valeurs > 1 correspondent à des hauteurs localisées sous 1,30 m.} & A & B & C \\
			\hline 
			$]0; 0.2[$ & 0.80(0.12) & 0.81(0.12) & 0.86(0.09) \\
			$[0.2; 0.4[$ & 0.68(0.13) & 0.67(0.13) & 0.71(0.12) \\
			$[0.4; 0.6[$ & 0.62(0.12) & 0.60(0.13) & 0.65(0.11) \\ 
			$[0.6; 0.8[$ & 0.62(0.12) & 0.61(0.14) & 0.65(0.14) \\
			$[0.8; 1[$ & 0.70(0.26) & 0.85(0.24) & 0.80(0.25) \\ 
			$> 1$ & 2.41(1.48) & 2.68(2.11) & 2.78(1.59) \\
			\hline
		\end{tabular}
		\end{minipage}
	\label{table_theta}
\end{table}

\subsection{Modélisation directe de la forme des tiges}

Dans cette partie c'est le coefficient de forme $\theta$ qui est modélisé avant d'être réinjecté dans l'équation (\ref{equa2}). On cherche à modéliser $\theta$ en fonction de la hauteur relative, au choix $x(h)=\frac{H-h}{H-1,3}$ ou $\frac{h}{H}$. La relation est plus simple à modéliser en fonction de $\frac{h}{H}$ car elle présente une forme relativement classique.

\begin{figure}[!h]
	\centering	
	\includegraphics[width=9cm]{theta_vs_h_rel_global.png}
	\caption{$\theta$ \textit{vs.} $h/H$ (une couleur par arbre).}
	\label{theta_vs_h_rel_global}
\end{figure}

Le modèle de base retenu est le suivant :
\begin{equation}
		\theta(\frac{h}{H}) = \exp(a-\frac{h}{H})^{b}+c\times(\frac{h}{H})^{d}+e
		\label{expo_variable_1}
\end{equation}
On fait de plus l'hypothèse que la forme du profil tend vers quelque chose de conique quand $h \to H$ (\textit{i.e.}, $\frac{h}{H} \to 1$) (en accord avec la littérature et de plus si on fait une petite imprécision dans cette zone ce n'est pas très grave car d'une part le volume représenté est faible et d'autre part cette zone n'est pas commercialisée) ce qui se traduit par $\theta(1)=1$ et qui permet de réduire le nombre de paramètres du modèle :
\begin{equation}
	\theta(\frac{h}{H}) = \exp(a-\frac{h}{H})^{b}+c\times(\frac{h}{H})^{d}+(1-\exp(a-1)^{b}-c)
	\label{expo_variable_2}
\end{equation}
où $a = 0.03234$, $b = 61.12$, $c = 0.3141$ et $d = 9.102$ (ajustement global \textit{via} la fonction \textit{nls} de R ; $n = 12878$ observations). 
\begin{figure}[!h]
	\centering	
	\includegraphics[width=9cm]{correlations_parametres.png}
	\caption{Corrélation entre les paramètres et l'âge et le $H/D$.}
	\label{correlations_parametres}
\end{figure}

\textit{En cours d'analyse et de rédaction...}


\subsection{Modélisation de la relation $AIB/AIC$ en attendant mieux... c'est-à-dire de travailler sur les profils complets}

En raison du besoin urgent dans le simulateur de croissance \textit{SimCop} d'une allocation annuelle plus juste de la matière ligneuse le long de la tige, nous avons travaillé sur la modélisation du ratio $AIB/AIC$, où $AIB$ est la surface du cerne mesurée à 1,30 m et $AIC$ est la surface du même cerne mesurée à la base du houppier. En effet, actuellement la surface du cerne est considérée constante sous le houppier (loi de Pressler) ce qui tend probablement à surestimer les $C_{130}$ des arbres dominés ou stressés.

Noël Le Goff et Jean-Marc Ottorini ont développé pour le Hêtre des modèles mettant en relation le ratio $AIB/AIC$ avec la longueur relative du houppier $LHP/HT$ (Equation \ref{AIB_AIC_Hetre_1} et \ref{AIB_AIC_Hetre_2}).

\begin{equation} 
	AIB/AIC = 1,45 \times (1 - \exp(-30 \times(LHP/HT)^3))
	\label{AIB_AIC_Hetre_1}
\end{equation}

\begin{equation} 
	AIB/AIC = 0,5 + 0,8 \times \exp(-10000 \times \exp(-30 \times (LHP/HT)))
	\label{AIB_AIC_Hetre_2}
\end{equation}

Nous avons cherché à retrouver le même type de relation chez les Douglas des jeux de données d'Ecouves puis d'Amance. Pour l'année 2010 à Ecouves et pour l'année 1991 à Amance, nous disposons de mesures de terrain pour la hauteur de la base du houppier prise au moment de l'abattage des arbres. Pour les autres années, nous avons reconstruit rétrospectivement la hauteur de la base du houppier à partir des valeurs de $Hf$ (voir section \ref{Modele_Deleuze}).

La Fig. \ref{AIB_AIC_vs_Longueur_relative_houppier} montre la relation entre $AIB/AIC$ et la hauteur relative de la base du houppier. Aucune tendance ne peut être observée. Le ratio $AIB/AIC$ n'apparaît corrélé à aucune des variables testées (Fig. \ref{pairs_AIB_AIC_vs_var_arbre}). Les Fig. \ref{AIB_AIC_vs_Annee_par_statut} et \ref{AIB_AIC_vs_Annee_par_traitement} montrent l'évolution du ratio $AIB/AIC$ avec le temps en fonction du statut social et du traitement sylvicole, respectivement, pour le jeu d'Ecouves. La Fig. \ref{AIB_AIC_vs_Age_par_site} montre l'évolution du ratio en fonction de l'âge des arbres pour les deux sites Ecouves (13 à 63 ans) et Amance (13 à 39 ans). 

En conclusion, nous pensons que dans notre cas il n'est pas pertinent de travailler sur ce ratio $AIB/AIC$ car à 1,30 m et à la base du houppier le profil de surface de cerne est très variable justement (Annexe \ref{Appendix_Shx}) et par conséquent ce ratio ne reflète absolument pas la partie du profil pour laquelle cette surface semble être plus ou moins constante et qui se situe bien au-dessus de 1,30 m. 



\begin{figure}[!h]
	\centering
	\includegraphics[width=\textwidth]{AIB_AIC_vs_Longueur_relative_houppier_comparaison_Ecouves_Amance.pdf}	
	\caption{$AIB/AIC$ en fonction de la longueur relative du houppier : en noir les estimations à partir des reconstructions rétrospectives de la remontée de la base du houppier ; en rouge les mesures lors de la dernière année au moment de l'abattage (2010 pour Ecouves et 1991 pour Amance).}
	\label{AIB_AIC_vs_Longueur_relative_houppier}
\end{figure}

\begin{figure}[!h]
	\centering
	%\includegraphics[width=\textwidth]{pairs_AIB_AIC_vs_var_arbre.pdf}	
	\includegraphics[width=\textwidth]{pairs_AIB_AIC_vs_var_arbre.png}	
	\caption{Corrélations entre les différentes variables niveau "arbre" et le ratio $AIB/AIC$ sur le jeu d'Ecouves uniquement.}
	\label{pairs_AIB_AIC_vs_var_arbre}
\end{figure}

\begin{figure}[!h]
	\centering
	\includegraphics[width=9cm]{AIB_AIC_vs_Annee_par_statut.pdf}	
	\caption{$AIB/AIC$ en fonction de l'année par statut social sur le jeu d'Ecouves uniquement.}
	\label{AIB_AIC_vs_Annee_par_statut}
\end{figure}

\begin{figure}[!h]
	\centering
	\includegraphics[width=9cm, page=1]{AIB_AIC_vs_Annee_par_traitement.pdf}	
	\caption{$AIB/AIC$ en fonction de l'année par traitement sylvicole sur le jeu d'Ecouves uniquement.}
	\label{AIB_AIC_vs_Annee_par_traitement}
\end{figure}

\begin{figure}[!h]
	\centering
	\includegraphics[width=9cm]{AIB_AIC_vs_Age_par_site.pdf}	
	\caption{$AIB/AIC$ en fonction de l'âge de l'arbre pour les deux jeux de données Ecouves et Amance.}
	\label{AIB_AIC_vs_Age_par_site}
\end{figure}

\FloatBarrier


\subsection{Modélisation des incréments annuels de la surface du cerne et confrontation à la loi de Pressler}

\label{Modele_Deleuze}

\subsubsection{Ajustement du modèle décrit par \cite{Deleuze2002}}

La fonction [$ShX$] retenue (équation \ref{ShX}) a été développée sur la base d'un modèle à base écophysiologique d'allocation du carbone. La fonction comporte 3 paramètres : $P_{r}$, $P_{f}$ et $z$. 

\begin{equation}
\Delta G(x) = x \frac{P_r \sinh(zx) + x P_f \sinh(z(H-x))}{\sinh(zH)}
\label{ShX}
\end{equation}
avec $\Delta G(x)$ l'incrément ligneux annuel en surface en m$^2$ à la position $x$ le long de la tige, $x$ étant la distance à l'apex en m et $H$ la hauteur totale de l'arbre en m.
Les dérivées première, seconde et troisième de cette fonction sont utilisées pour calculer des indicateurs $H_f$, $H_s$ et $H_b$ de la hauteur de la base du houppier fonctionnel selon le calcul détaillé en Annexe \ref{Appendix_Shx}. Nous disposons de mesures de terrain de hauteur de la base du houppier, de hauteur de la 1\iere branche verte et de hauteur du 1\ier frottis uniquement pour la dernière année avant abattage (2010 pour Ecouves et 1991 pour Amance). La figure \ref{ShX_et_derivees_par_arbre} montre quelques exemples de résultats obtenus (voir l'ensemble des résultats en Annexe \ref{Appendix_Shx}). L'Annexe \ref{Appendix_Shx_2} présente les ajustements du modèle sur les profils d'incréments des années passées : 1960, 1970, 1980, 1990, 2000 et 2010 pour Ecouves et 1965, 1970, 1975, 1980, 1985 et 1990 pour Amance. Pour le site d'Ecouves, le modèle semble presque mieux se comporter pour les années passées qu'en 2010.  

\begin{figure}[!h]
	\centering
	\includegraphics[scale=.4, page=1339]{ShX_et_derivees_par_arbre_par_annee_Ecouves}	
	\includegraphics[scale=.4, page=1354]{ShX_et_derivees_par_arbre_par_annee_Ecouves}	
	\includegraphics[scale=.4, page=1387]{ShX_et_derivees_par_arbre_par_annee_Ecouves}	
	\caption{Dérivées première, seconde et troisième de la fonction [$ShX$] obtenues pour trois années de l'arbre 27 d'Ecouves (co-dominant, densité faible).}
	\label{ShX_et_derivees_par_arbre}
\end{figure}

\FloatBarrier
\bigskip

Nous avons analysés les résidus des ajustements du modèle de profil de surface de cerne. Aucune différence significative sur la moyenne des résidus des surfaces de cernes n'a pu être mise en évidence entre les deux sites. 

Sur le site d'Ecouves, aucune différence n'a pu être observée entre les différents statuts sociaux, ni entre les différents traitements A, B et C.

Sur le site d'Amance, le traitement D4 apparaît significativement différent des autres (B1, C3 et D1) en terme de moyenne de résidus.

La Fig. \ref{residus_modele_ini_vs_hauteur} illustre l'évolution des résidus avec la hauteur dans l'arbre. Pour le site d'Amance aucun problème particulier ne semble à signaler. En revanche pour le site d'Ecouves nous observons une surestimation de la surface de cerne par le modèle à l'empattement et une sous-estimation dans le houppier. Le modèle bien que relativement souple ne permet pas de bien coller aux formes de profils rencontrées à Ecouves. La surestimation à l'empattement est déjà observée avant 39 ans sur les arbres d'Ecouves donc il s'agit probablement d'un effet du site plutôt que de la différence d'âge entre les deux sites. En revanche, la sous-estimation des surfaces de cernes dans le houppier ne semble se produire que pour les arbres plus âgés d'Ecouves. 

\begin{figure}[!h]
	\centering
	\includegraphics[width=\textwidth]{residus_ajustement_initial_vs_hauteur.pdf}	
	\caption{Représentation des résidus en fonction de la hauteur dans l'arbre pour chacun des deux sites. La hauteur 1,30 m est indiquée par une ligne en pointillés grise. Les moyennes des résidus pour les hauteurs communes de 30 cm, 80 cm et 130 cm pour Ecouves et de 130 cm pour Amance sont indiquées par des points rouges. Trente-cinq résidus supérieurs en valeur absolue à 0.002 m$^{2}$ ne sont pas représentés pour le site d'Ecouves.}
	\label{residus_modele_ini_vs_hauteur}
\end{figure}

\FloatBarrier
\bigskip 

La figure \ref{correlations_base_houppier} montre qu'il existe une bonne relation entre la base du houppier mesurée sur le terrain au moment de l'abattage des arbres et les indicateurs $H_f$ et $H_s$, notamment.

$H_f$ a le mérite d'être moins sensible que $H_s$ et par conséquent nous l'avons retenu pour la prédiction de la hauteur de base du houppier (Équation \ref{pred_base_houppier}). Plus exactement c'est une valeur de $H_f$ lissé qui est utilisée dans l'équation de façon à s'assurer d'une remontée en cas de problème une année donnée.

\begin{multline} 
	Hauteur\_base\_houppier(m) = -5,17446 + 0,78679 \times H_f\_smoothed(m) + \\ 0,04652 \times H/D
	\label{pred_base_houppier}
\end{multline} 

\begin{figure}[!h]
	\centering
	\includegraphics[width=\textwidth]{Correlations.pdf}	
	\caption{Corrélations entre la base du houppier mesurée sur le terrain (en 2010 pour le site d'Ecouves et en 1991 pour le site d'Amance) et les indicateurs obtenus avec la fonction [$ShX$].}
	\label{correlations_base_houppier}
\end{figure}

La figure \ref{Remontees_H_Hf_Hs_Hb} montre les remontées théoriques de la base du houppier obtenues avec les 3 indicateurs calculés sur les profils de cernes de chaque année pour le même arbre. Les résultats obtenus pour l'ensemble des arbres sont présentés en Annexe \ref{Appendix_Remontees}.

\begin{figure}[!h]
	\centering
	\includegraphics[width=9cm]{Exemple_remontees_H_Hf_Hs_Hb.pdf}
	\caption{Remontée théorique de la base du houppier obtenue avec la fonction [$ShX$] pour l'arbre 27 d'Ecouves.}
	\label{Remontees_H_Hf_Hs_Hb}
\end{figure}

\begin{figure}[!h]
	\centering
	\includegraphics[width=9cm]{Hbh_predite_vs_Htot_par_site.pdf}
	\caption{Relation entre hauteurs totales et hauteurs de base de houppier prédites et mesurées (sites d'Ecouves et d'Amance).}
	\label{Hbh_predite_vs_Htot_par_site}
\end{figure}
\FloatBarrier

\subsubsection{Modélisation des profils d'incréments annuels à partir de variables "arbre"}

La Fig. \ref{pairs_Pr_Pf_z_vs_var_arbre} montre les corrélations existant entre les 3 paramètres du modèle (ajustement pour l'année $n$) et des variables intrinsèques à l'arbre (pour l'année $n-1$). 

\begin{figure}[!h]
	\centering
	%\includegraphics[width=\textwidth]{pairs_Pr_Pf_z_vs_var_arbre.pdf}
	\includegraphics[width=\textwidth]{pairs_Pr_Pf_z_vs_var_arbre.png}
	\caption{Corrélations entre les différentes variables niveau "arbre" et les paramètres $P_{r}$, $P_{f}$ et $z$ du modèle (sites d'Ecouves et d'Amance).}
	\label{pairs_Pr_Pf_z_vs_var_arbre}
\end{figure}

Nous avons choisi de modéliser $z$ en fonction du $D_{130}$ de l'année $n-1$ (Équation \ref{equa_modelisation_z}). La fonction décroit vers une asymptote $a$ quand le diamètre tend vers l'infini ce qui assure que $z$ reste positif $\forall D130$.

\begin{equation} 
	z = a+\exp(b-D130_{n-1})^c
	\label{equa_modelisation_z}
\end{equation}

avec $a = 0,08123$, $b = -17,07634$ et $c = 0,05115$.

\begin{figure}[!h]
	\centering	
	\includegraphics[width=9cm]{modelisation_z.pdf}
	\caption{Modélisation du paramètre $z$ en fonction du $D_{130}$ de l'année $n-1$.}
	\label{modelisation_z}
\end{figure}

\begin{figure}[!h]
	\centering	
	\includegraphics[width=\textwidth]{modelisation_z_residus.pdf}
	\caption{Résidus du modèle en fonction de variables "arbre".}
	\label{modelisation_z_residus}
\end{figure}

\FloatBarrier

Nous avons choisi de modéliser $P_{f}$ en fonction de l'âge de l'année $n-1$ (modèle 1, équation \ref{equa_modelisation_Pf_1}). Nous avons également testé un modèle à deux variables avec le $D_{130}$ de l'année $n-1$ (modèle 2, équation \ref{equa_modelisation_Pf_2}).

\begin{equation} 
	P_{f} = a+\exp(b-Age_{n-1})^c
	\label{equa_modelisation_Pf_1}
\end{equation}

avec $a = -7,139e-05$, $b = -2,888e+02$ et $c = 2,611e-02$. 

\begin{equation} 
	P_{f} = a+\exp(b-Age_{n-1})^c + d \times D130_{n-1} 
	\label{equa_modelisation_Pf_2}
\end{equation}

avec $a = -1,774e-04$, $b = -2,676e+02$, $c = 2,739e-02$ et $d = 2,629e-06$. 

La correction avec le $D_{130}$ permet de corriger légèrement les tendances observées dans les résidus. L'amélioration n'est cependant pas suffisante pour justifier le paramètre supplémentaire. Nous utiliserons donc par la suite le modèle 1 (équations \ref{equa_modelisation_Pf_1}).

\begin{figure}[!h]
	\centering	
	\includegraphics[width=9cm]{modelisation_Pf.pdf}
	\caption{Modélisation du paramètre $P_{f}$ en fonction de l'âge de l'arbre l'année $n-1$.}
	\label{modelisation_Pf}
\end{figure}

\begin{figure}[!h]
	\centering	
	\includegraphics[width=\textwidth]{modelisation_Pf_residus_modele_1.pdf}
	\caption{Résidus du modèle 1 (équation \ref{equa_modelisation_Pf_1}) en fonction de variables "arbre".}
	\label{modelisation_Pf_residus_1}
\end{figure}

\begin{figure}[!h]
	\centering	
	\includegraphics[width=\textwidth]{modelisation_Pf_residus_modele_2.pdf}
	\caption{Résidus du modèle 2 (équation \ref{equa_modelisation_Pf_2}) en fonction de variables "arbre".}
	\label{modelisation_Pf_residus_2}
\end{figure}

\FloatBarrier

Enfin nous avons choisi de modéliser $P_{r}$ en fonction de la largeur moyenne à 1,30 m des 5 derniers cernes pour l'année $n-1$ (modèle 1, équation \ref{equa_modelisation_Pr_1}).

\begin{equation} 
	P_{r} = a \times LC5_{n-1}
	\label{equa_modelisation_Pr_1}
\end{equation} 

avec $a = 3,746e-05$.

Avec cette formulation $P_{r}$ est toujours $> 0$. 
Nous avons également testé un modèle à trois paramètres avec une constante et le $D_{130}$ de l'année $n-1$ (modèle 2, équation \ref{equa_modelisation_Pr_2}).

\begin{equation} 
	P_{r} = a \times LC5_{n-1} + b + c \times D130_{n-1} 
	\label{equa_modelisation_Pr_2}
\end{equation} 

avec $a = 3,820e-05$, $b=-2,201e-05$ et $c=8,729e-07$.

\begin{figure}[!h]
	\centering	
	\includegraphics[width=9cm]{modelisation_Pr.pdf}
	\caption{Modélisation du paramètre $P_{r}$ en fonction de la largeur moyenne à 1,30 m des 5 derniers cernes à l'année $n-1$.}
	\label{modelisation_Pr}
\end{figure}

\begin{figure}[!h]
	\centering	
	\includegraphics[width=\textwidth]{modelisation_Pr_residus_modele_1.pdf}
	\caption{Résidus du modèle 1 (équation \ref{equa_modelisation_Pr_1}) en fonction de variables "arbre".}
	\label{modelisation_Pr_residus_1}
\end{figure}

\begin{figure}[!h]
	\centering	
	\includegraphics[width=\textwidth]{modelisation_Pr_residus_modele_2.pdf}
	\caption{Résidus du modèle 2 (équation \ref{equa_modelisation_Pr_2}) en fonction de variables "arbre".}
	\label{modelisation_Pr_residus_2}
\end{figure}

%VO : Pour le site d'Amance il semble y avoir une année présentant de fort résidus pour des valeurs de $LC5_{n-1}$ a priori importante pour quelques années particulières.
La forme des nuages obtenus dans les figures \ref{modelisation_Pr_residus_1} et \ref{modelisation_Pr_residus_2} montre une hétéroscédasticité des résidus, notamment avec la largeur des cernes et le rapport $H/D$.
En particulier, on observe de forts résidus sur le site d'Amance pour quelques années de forte croissance, lorsque les arbres avaient entre 10 et 20 ans.

La prise en compte du $D_{130}$ n'améliore pas sensiblement les résultats. Nous utiliserons donc par la suite le modèle 1 (équation \ref{equa_modelisation_Pr_1}).

\FloatBarrier

\bigskip

Nous avons ensuite ajusté le modèle global en recherchant un effet site sur les paramètres du modèle avec la fonction \textit{nlsList}. La Fig. \ref{nlsList_effet_site} montre qu'il y a un effet site sur les 7 paramètres du modèle.

\begin{equation}
\Delta G(x) (m^{2}) = x \frac{P_r \sinh(zx) + x P_f \sinh(z(H-x))}{\sinh(zH)}
\label{modelisation_surface}
\end{equation}

où $x$ (m) est la distance à l'apex et $H$ (m) est la hauteur totale de l'arbre

et avec : 

\begin{equation}
	z = a+\exp(b-D130_{n-1})^c
\end{equation}

\begin{equation}
	P_{f} = d+\exp(e-Age_{n-1})^f
\end{equation}

\begin{equation}
	P_{r} = g \times LC5_{n-1}
\end{equation}

\begin{figure}[!h]
	\centering	
	\includegraphics[width=\textwidth]{nlsList_effet_site.pdf}
	\caption{Recherche d'un effet site sur les paramètres du modèle.}
	\label{nlsList_effet_site}
\end{figure}

Plutôt que de rentrer l'effet site en effet fixe dans le modèle nous avons choisi de le rentrer en effet aléatoire sur chacun des paramètres de façon à obtenir le meilleur modèle moyen et à pouvoir appliquer, le cas échéant, la partie à effets fixes du modèle sur un nouveau site. 

Après ajustement, le paramètre $d$ devient non significatif. Nous décidons cependant de le garder dans le modèle.
% Fred : il faudrait expliquer pourquoi !!

\FloatBarrier

La Fig. \ref{modelisation_globale_residus} montre les résidus du modèle en fonction de différentes variables. Nous pouvons observer une dérive surtout visible pour les grandes surfaces de cernes. Ces forts résidus sont obtenus pour des hauteurs relatives faibles, donc à l'empattement, et pour des arbres plutôt âgés du site d'Ecouves. Cela signifie que le modèle ne s'ajuste pas parfaitement aux données, notamment à l'empattement et qu'il sera probablement difficile d'estimer correctement un diamètre à 1,30 m à partir de ces prédictions. Toutefois cette dérive est à relativiser car le nombre de surfaces de cernes observées au-delà de 5000 mm$^{2}$ est complètement négligeable (Fig. \ref{effectifs_surfaces_cernes}). 

\begin{figure}[!h]
	\centering	
	\includegraphics[width=\textwidth]{modelisation_globale_residus_effet_aleatoire_site.png}
	\caption{Résidus du modèle global (Équation \ref{modelisation_surface}) avec des effets aléatoires site sur les paramètres du modèle.}
	\label{modelisation_globale_residus}
\end{figure}

\begin{figure}[!h]
	\centering	
	\includegraphics[width=9cm]{Effectifs_par_surface_cerne.pdf}
	\caption{Effectifs en fonction des surfaces de cernes.}
	\label{effectifs_surfaces_cernes}
\end{figure}

Les Fig. \ref{surf130_obs_vs_surf130_pred} et \ref{volume_obs_vs_volume_pred} représentent les valeurs observées en fonction des valeurs prédites pour la surface de cerne à 1,30 m et pour le volume total du cerne. La RMSE pour la surface de cerne à 1,30 m est de $6,45 \times 10^{-4}$ m$^{2}$ et la RMSE relative correspondante est de 30 \%. 
La RMSE à Ecouves est de $5,83 \times 10^{-4}$ m$^{2}$ tandis qu'à Amance elle est de $7,48 \times 10^{-4}$ m$^{2}$. 

L'erreur moyenne est de $-2,60 \times 10^{-4}$ m$^{2}$ à Ecouves tandis qu'elle est de $+4,96 \times 10^{-6}$ m$^{2}$ à Amance. Il y a donc une surestimation de la surface de cerne à 1,30 m à Ecouves et une sous-estimation mais de moindre importance à Amance.

La RMSE pour le volume total du cerne est de $7,58 \times 10^{-3}$ m$^{3}$ et la RMSE relative correspondante est de 23 \%.

\begin{figure}[!h]
	\centering	
	\includegraphics[width=\textwidth]{surf130_obs_vs_surf130_pred.pdf}
	\caption{Surface de cerne observée à 1,30 m en fonction de la surface de cerne prédite par le modèle.}
	\label{surf130_obs_vs_surf130_pred}
\end{figure}

\begin{figure}[!h]
	\centering	
	\includegraphics[width=\textwidth]{volume_obs_vs_volume_pred.pdf}
	\caption{Volume du cerne observé en fonction du volume du cerne prédit par le modèle.}
	\label{volume_obs_vs_volume_pred}
\end{figure}

La Fig. \ref{modelisation_croissance_arbre_moyen_1} présente l'évolution des résidus pour l'estimation du $D_{130}$ au cours du temps lorsqu'on modélise l'arbre moyen. Le calcul est fait en cumulant au cours du temps les diamètres à 1,30 m modélisés. Ici le modèle prend en entrée le vrai diamètre à 1,30 m de l'année précédente ($D130_{n-1}$) et la vraie largeur moyenne à 1,30 m des 5 derniers cernes de l'année précédente ($LC5_{n-1}$). Au travers de cette dernière variable le modèle prend partiellement en compte le climat et l'effet des éclaircies. Les résultats pour Ecouves montrent une tendance à l'augmentation de la surestimation du modèle avec l'âge. 

La Fig. \ref{modelisation_croissance_arbre_moyen_2} présente les résidus que l'on obtiendrait si on utilisait en entrée du modèle les sorties calculées à partir des années précédentes. Les résidus sont beaucoup trop importants. La dérive pour Ecouves est très forte. Le point positif est que pour les arbres d'Amance, bien que les résidus soient importants, il n'y a pas de dérive marquée. 

\begin{figure}[!h]
	\centering	
	\includegraphics[width=\textwidth]{modelisation_croissance_arbre_moyen_1.pdf}
	\caption{Résidus de la croissance en $D_{130}$ de l'arbre moyen modélisé à partir des variables mesurées.}
	\label{modelisation_croissance_arbre_moyen_1}
\end{figure}

\begin{figure}[!h]
	\centering	
	\includegraphics[width=\textwidth]{modelisation_croissance_arbre_moyen_2.pdf}
	\caption{Résidus de la croissance en $D_{130}$ de l'arbre moyen modélisé à partir de variables prédites (processus itératif).}
	\label{modelisation_croissance_arbre_moyen_2}
\end{figure}

%Tout dépend de la façon dont on souhaiterait appliquer un tel modèle. Il serait en effet possible de mettre en entrée du modèle la hauteur totale de l'arbre plutôt que le $D_{130}$ ce qui éviterait d'avoir à utiliser notre modèle, relativement peu précis à 1,30 m, de façon itérative. 
La dérive observée s'explique par l'accumulation des erreurs sur le $D_{130}$ qui est estimé par le modèle (avec une erreur relativement forte) et sert lui même d'entrée pour l'année suivante.
Nous présentons dans le paragraphe suivant un modèle alternatif permettant de remédier à ce problème.

% Fred : comme on vient de dire que le modèle n'était pas bon, c'est un peu bête de continuer avec...
% 	Je propose de virer la suite jusqu'au § "Recherche d'un modèle équivalent sans dérive"
% 	et de refaire les analyses des résidus / année avec le modèle sans dérive

La Fig. \ref{modelisation_globale_residus_vs_annee} montre les résidus du modèle en fonction de l'année calendaire. Il est intéressant de constater que le modèle surestime la surface du cerne pour les années de sécheresse. Effectivement pour ces années caractéristiques les cernes sont réduits mais le modèle n'ayant pas d'entrées climatiques correctes ne peut pas en tenir compte. Le $LC5_{n-1}$ est contient toutefois une information climatique mais qui intervient avec un décalage dans le temps.   

De la même manière, le modèle dans sa forme actuelle ne peut pas modéliser la réaction de l'arbre après une éclaircie. Dans ce cas le modèle sous-estimera la croissance de l'arbre. Cependant comme le modèle prend en entrée la largeur moyenne des 5 derniers cernes à 1,30 m (de l'année précédente) il est possible de prendre une partie de l'effet en compte. Mais le modèle sera incapable de modéliser correctement la 1ère année après l'éclaircie de cette façon. 

Au final on pourra donc avoir une surestimation de petites surfaces de cernes (années climatiques particulière ou forte compétition) et une sous-estimation de grandes surfaces de cernes (après les éclaircies notamment) et cela peut contribuer également à la tendance observée dans les résidus en fonction de la surface de cerne.

\bigskip 

\textbf{N.B. Comme perspectives à court terme nous avons de réajuster un modèle moyen de surfaces de cernes ne prenant pas du tout en compte le climat, ni l'effet des éclaircies, donc sans $LC5_{n-1}$, uniquement âge et taille (voire forme) de l'arbre, ces dernières étant des variables obtenues par cumul et donc non sensibles aux perturbations de hautes fréquences, afin de mieux faire ressortir l'effet d'années caractéristiques sur les profiles des résidus. Les données climatiques sont en cours d'acquisition.}  

\bigskip

\begin{figure}[!h]
	\centering	
	\includegraphics[width=\textwidth]{modelisation_globale_residus_effet_aleatoire_site_vs_annee.pdf}
	\caption{Résidus (moyenne et intervalle de confiance) du modèle global (Équation \ref{modelisation_surface}) avec des effets aléatoires "site" sur les paramètres du modèle en fonction de l'année calendaire. Les lignes verticales grises en pointillés correspondent aux dates des éclaircies. Les lignes verticales rouges en pointillés correspondent à des années de sécheresse.}
	\label{modelisation_globale_residus_vs_annee}
\end{figure}

Par ailleurs, pour l'année de sécheresse de 1976 nous avons tracé les résidus en fonction de la hauteur relative dans l'arbre (Fig. \ref{modelisation_globale_residus_vs_hauteur_relative_1976}) pour vérifier la concordance avec des observations de la littérature qui relatent qu'en cas de stress intense les profils d'incréments annuels sont plus réduit dans le bas des tiges (avec même parfois pour certaines essences comme le Hêtre la disparition complète du cerne) que dans le haut (allocation prioritaire en haut).

\begin{figure}[!h]
	\centering	
	\includegraphics[width=\textwidth]{residus_vs_hauteur_relative_1976.pdf}
	\caption{Résidus en fonction de la hauteur relative pour l'année de sécheresse de 1976.}
	\label{modelisation_globale_residus_vs_hauteur_relative_1976}
\end{figure}

Dans la Fig. \ref{residus_vs_annee_par_hauteur}, nous avons représenté les résidus par traitement sylvicole (A = témoin, B et C) et par hauteur dans l'arbre (0,3 m, 1,30 m et la rondelle située sous la base du houppier quand celle-ci était différente de celle à 1,30 m). 
On observe une plus grande variabilité des résidus à l'empattement ce qui n'est pas très étonnant étant donné que c'est dans cette partie que la surface du cerne est la plus importante en même temps qu'elle est très variable. 
Les surestimations du modèle pour 1976 et au début des années 1990 (sécheresses recensées mais à confirmer avec des données climatiques locales) sont présentent dans tous les cas de figure. La surestimation est d'autant plus forte qu'on se situe vers le bas de la tige ce qui confirme le changement d'allocation lié à un stress ressenti ces années là avec une allocation plus réduite vers le bas que vers le haut de la tige. 
Tous les peuplements présentent un schéma similaire à partir de l'année 1995 (décroissance de la courbe des résidus), ce qui semble être lié à un évènement climatique à élucider.
Les résidus sont moins variables dans le peuplement témoin A et la variabilité additionnelle observée dans les peuplements éclaircis serait à attribuer aux interventions sylvicoles.
Une reprise de croissance qui aurait lieu dans tous les traitements pourrait être attribuer à une tempête. Une reprise de croissance qui n'aurait lieu que dans les peuplements B et/ou C serait à attribuer aux éclaircies.  

\begin{figure}[!h]
	\centering	
	\includegraphics[width=\textwidth]{residus_vs_annee_par_hauteur.pdf}
	\caption{Résidus (moyenne et intervalle de confiance) du modèle global (Équation \ref{modelisation_surface}) avec des effets aléatoires "site" sur les paramètres du modèle en fonction de l'année calendaire. Les lignes verticales grises en pointillés correspondent aux dates des éclaircies. Les lignes verticales rouges en pointillés correspondent à des années de sécheresse (à vérifier). Les données sont présentées par traitement sylvicole et pour trois hauteurs : 30 cm, 1,30 m et la rondelle située juste sous la base du houppier estimée. Seul le site d'Ecouves est représenté car sur le site d'Amance il n'y a pas de rondelles prélevées sous 1,30 m.}
	\label{residus_vs_annee_par_hauteur}
\end{figure}

Les résidus par site, par arbre et par année sont présentés en Annexe \ref{Residus_modele_increment_par_arbre_par_annee}. Les arbres présentant les plus forts résidus sont deux arbres dominants du traitement intermédiaire (B) du site d'Ecouves (arbres 1 et 4).

\FloatBarrier

\subsubsection{Recherche d'un modèle équivalent sans dérive...}

Dans la section précédente nous avions choisi les meilleures variables prédictives pour chacun des paramètres $z$, $P_f$ et $P_r$. Parmi elles figuraient le $D130_{n-1}$ et la largeur moyenne des 5 derniers cernes à 1,30 m ($LC5_{n-1}$). Pour l'application du modèle à partir de variables mesurées sur le terrain ou estimées correctement par ailleurs il n'y a pas de problème mais si on itère et qu'on ré-applique notre modèle en utilisant les sorties de l'année précédente pour estimer $D130_{n-1}$ et $LC5_{n-1}$ alors une dérive est observée car le modèle est fait pour prédire des volumes de cernes et non des mesures à 1,30 m, zone où de plus la surface de cerne est très variable.  

Pour cela nous avons cherché à calculer des variables équivalentes à partir des volumes. Dans la liste des variables candidates, nous avons ajouté : le volume du cerne de l'année précédente ($V\_cerne_{n-1}$), le volume du tronc de l'année précédente, un pseudo-équivalent du $D130_{n-1}$ calculé comme étant  $D\_volume_{n-1} = 2 \times \sqrt{1/\pi \times V\_cerne_{n-1}/H_{n-1}}$, un pseudo-équivalent du $H/D_{n-1}$ calculé comme étant $H/D\_volume_{n-1} = H_{n-1}/D\_volume_{n-1}$ et un pseudo-équivalent de la largeur moyenne des 5 derniers cernes à 1,30 m ($LC5\_volume_{n-1}$) calculée à partir de soustractions des pseudo $D_{130}$. Il existe des relations linéaires très fortes entre les différentes variables :
\label{definitions_variables_volume}

\begin{equation}
D130 = -1,19048 + 1,53221 \times D\_volume    
\end{equation}

avec un R$^{2}$ de 0,99.

\begin{equation}
H/D = 7,81000 + 0,64418 \times H/D\_volume    
\end{equation}

avec un R$^{2}$ de 0,85.

\begin{equation}
LC5 = 0,28185 + 0,70962 \times LC5\_volume    
\end{equation}

avec un R$^{2}$ de 0,91.

\bigskip

La Fig. \ref{pairs_Pr_Pf_z_vs_var_arbre_2} montre les corrélations existant entre les 3 paramètres du modèle (ajustement pour l'année $n$) et des variables intrinsèques à l'arbre (pour l'année $n-1$) incluant les nouvelles variables calculées. Comme déjà montré ci-dessus, nos nouvelles variables $D\_volume_{n-1}$, $H/D\_volume_{n-1}$ et $LC5\_volume_{n-1}$ sont très corrélées avec celles qu'elles sont censées remplacer. 

\begin{figure}[!h]
	\centering
	\includegraphics[width=\textwidth]{pairs_Pr_Pf_z_vs_var_arbre_2.png}
	\caption{Corrélations entre les différentes variables niveau "arbre" et les paramètres $P_{r}$, $P_{f}$ et $z$ du modèle (sites d'Ecouves et d'Amance).}
	\label{pairs_Pr_Pf_z_vs_var_arbre_2}
\end{figure}

\bigskip 

Nous avons choisi de modéliser $z$ en fonction du $D\_volume_{n-1}$ de l'année $n-1$ estimé à partir du volume du tronc (Équation \ref{equa_modelisation_z_sans_derive}). La corrélation de Spearman entre le paramètre $z$ et $D\_volume_{n-1}$ est de -0,80. 

\begin{equation} 
	z = a+\exp(b-D\_volume_{n-1}(cm))^c
	\label{equa_modelisation_z_sans_derive}
\end{equation}

avec $a = 0,07603$, $b = -11,33489$ et $c = 0,07383$.

\begin{figure}[!h]
	\centering	
	\includegraphics[width=9cm]{modelisation_z_sans_derive.pdf}
	\caption{Modélisation du paramètre $z$ en fonction du $D_{130}$ de l'année $n-1$.}
	\label{modelisation_z_sans_derive}
\end{figure}

\begin{figure}[!h]
	\centering	
	\includegraphics[width=\textwidth]{modelisation_z_residus_sans_derive.pdf}
	\caption{Résidus du modèle en fonction de variables "arbre".}
	\label{modelisation_z_residus_sans_derive}
\end{figure}

\FloatBarrier

\bigskip

Nous avons choisi de modéliser $P_f$ en fonction de l'âge de l'année précédente (Équation \ref{equa_modelisation_Pf_1_sans_derive}). La corrélation de Spearman entre le paramètre $P_f$ et $Age_{n-1}$ est de -0,82.
% Fred : sauf que précédemment on n'avait pas gardé le $D_{130}$ :
Comme précédemment, faire intervenir le $D\_volume_{n-1}$ dans l'équation permet de corriger une tendance à la sous-estimation des surfaces de cernes pour les gros diamètres.

L'équation de type Gompertz négatif utilisée ici est plus adaptée que celle utilisée dans le modèle précédent (avec dérive) et on pourrait éventuellement améliorer le modèle développé précédemment.

\begin{equation} 
	P_{f} = a \times \exp(b \times \exp(c \times Age_{n-1}))+d \times D\_volume_{n-1}(cm) 
	\label{equa_modelisation_Pf_1_sans_derive}
\end{equation}

avec $a = 0.00046$, $b = -0.24279$, $c = 0.05450$ et $d = 2,676e-06$. 

\begin{figure}[!h]
	\centering	
	\includegraphics[width=9cm]{modelisation_Pf_sans_derive.pdf}
	\caption{Modélisation du paramètre $P_f$ en fonction de l'âge de l'année $n-1$.}
	\label{modelisation_P_f_sans_derive}
\end{figure}

\begin{figure}[!h]
	\centering	
	\includegraphics[width=\textwidth]{modelisation_Pf_residus_sans_derive.pdf}
	\caption{Résidus du modèle en fonction de variables "arbre".}
	\label{modelisation_P_f_residus_sans_derive}
\end{figure}

\FloatBarrier

\bigskip

Nous avons choisi de modéliser $P_r$ en fonction du $LC5\_volume_{n-1}$ avec également un terme correctif fonction du $D\_volume_{n-1}$ (Équation \ref{equa_modelisation_Pr_sans_derive}). La corrélation de Spearman entre le paramètre $P_r$ et la variable $LC5\_volume_{n-1}$ est de 0,81.  

\begin{equation} 
	P_{r} = a + b \times LC5\_volume_{n-1}(mm) + c \times D\_volume_{n-1}(cm)  
	\label{equa_modelisation_Pr_sans_derive}
\end{equation}

avec $a = -3.11337e-05$, $b = 3.02674e-05$ et $c = 1.69586e-06$.

\begin{figure}[!h]
	\centering	
	\includegraphics[width=9cm]{modelisation_Pr_sans_derive.pdf}
	\caption{Modélisation du paramètre $P_r$ en fonction de $LC5\_volume_{n-1}$.}
% Fred : on ne parle nulle part du modèle 1 sans d_vol_prec
	\label{modelisation_Pr_sans_derive}
\end{figure}

\begin{figure}[!h]
	\centering	
	\includegraphics[width=\textwidth]{modelisation_Pr_residus_sans_derive.pdf}
	\caption{Résidus du modèle en fonction de variables "arbre".}
	\label{modelisation_Pr_residus_sans_derive}
\end{figure}

\FloatBarrier

Nous avons ensuite ajusté le modèle global en recherchant un effet site (problèmes de convergence pour la recherche d'un effet "traitement" et "arbre") sur les paramètres du modèle avec la fonction \textit{nlsList}.

\begin{equation}
\Delta G(x) (m^{2}) = x \frac{P_r \sinh(zx) + x P_f \sinh(z(H-x))}{\sinh(zH)}
\label{modelisation_surface_sans_derive}
\end{equation}

où $x$ (m) est la distance à l'apex et $H$ (m) est la hauteur totale de l'arbre

et avec : 

\begin{equation}
	z = a+\exp(b-D\_volume_{n-1})^c
\end{equation}

\begin{equation}
	P_{f} = d \times \exp(e \times \exp(f \times Age_{n-1}))+g \times D\_volume_{n-1}(cm) 
\end{equation}

\begin{equation}
	P_{r} = h + i \times LC5\_volume_{n-1}(mm) + j \times D\_volume_{n-1}(cm)
\end{equation}
       
Dans un premier temps, pour des raisons de problèmes de convergence, nous allons fixer les paramètres $g$ et $j$ à 0 et $d$ (non significatif dans le modèle global) à sa valeur obtenue lors de l'ajustement local ($d = 0.0016331$) lorsque $g$ était pris égal à 0 (non présenté). La Fig. \ref{nlsList_effet_site_sans_derive} donne les résultats de la \textit{nlsList} pour la recherche d'un effet site sur les paramètres du modèle. L'effet site est significatif sur les paramètres $e$, $f$, $h$ et $i$.
         
\begin{figure}[!h]
	\centering	
	\includegraphics[width=\textwidth]{nlsList_effet_site_sans_derive.pdf}
	\caption{Recherche d'un effet site sur les paramètres du modèle.}
	\label{nlsList_effet_site_sans_derive}
\end{figure}

Plutôt que de rentrer l'effet site en effet fixe dans le modèle nous avons choisi de le rentrer en effet aléatoire sur les paramètres $e$, $f$, $h$ et $i$ concernés de façon à obtenir le meilleur modèle moyen et à pouvoir appliquer, le cas échéant, la partie à effets fixes du modèle sur un site autre que les deux utilisés pour ajuster le modèle. Dans le modèle mixte nous avons ensuite tenté de faire rentrer à nouveau les paramètres $g$ et $j$ qui se sont avérés être significatifs. Cependant, un modèle sans ces deux paramètres est également possible. 

Les valeurs estimées des paramètres (effets fixes du modèle décrit ci-dessus) sont les suivantes ($d$ étant par ailleurs fixé à 0,0016331) : 

$a  =  8,33718e-02$, $b = -7,17670e+00$, $c = 8,84324e-02$, $e = -1,05642e+00$, $f = 2,95686e-02$, $g = 2,01150e-06$, $h = -5,76793e-05$, $i=3,34321e-05$ et $j=2,36639e-06$. 

\begin{figure}[!h]
	\centering	
	%\includegraphics[width=\textwidth]{modelisation_globale_residus_effet_aleatoire_site_sans_derive.pdf}
	% (conversion par pdftoppm puis convert)
	\includegraphics[width=\textwidth]{modelisation_globale_residus_effet_aleatoire_site_sans_derive}
	\caption{Résidus du modèle global (Équation \ref{modelisation_surface_sans_derive}) avec des effets aléatoires site sur les paramètres $d$, $g$ et $h$ du modèle.}
	\label{modelisation_globale_residus_sans_derive}
\end{figure}

La Fig. \ref{modelisation_globale_residus_sans_derive} montre les résidus du modèle global avec effet aléatoire "Site" sur 4 paramètres. Les résultats obtenus semblent assez similaires à ceux précédemment obtenus.

Les Fig. \ref{surf130_obs_vs_surf130_pred_sans_derive} et \ref{volume_obs_vs_volume_pred_sans_derive} représentent les valeurs observées en fonction des valeurs prédites pour la surface de cerne à 1,30 m et pour le volume total du cerne. La RMSE pour la surface de cerne à 1,30 m est de $6,68 \times 10^{-4}$ m$^{2}$ et la RMSE relative correspondante est de 31 \%. La RMSE à Ecouves est de $6,20 \times 10^{-4}$ m$^{2}$ tandis qu'à Amance elle est de $7,52 \times 10^{-4}$ m$^{2}$. 

L'erreur moyenne est de $-1,93 \times 10^{-4}$ m$^{2}$ à Ecouves tandis qu'elle est de $2,76 \times 10^{-5}$ m$^{2}$ à Amance. Il y a donc une surestimation de la surface de cerne à 1,30 m à Ecouves et une sous-estimation mais de moindre importance à Amance.

La RMSE pour le volume total du cerne est de $7,22 \times 10^{-3}$ m$^{3}$ et la RMSE relative correspondante est de 22 \%.

\begin{figure}[!h]
	\centering	
	\includegraphics[width=\textwidth]{surf130_obs_vs_surf130_pred_sans_derive.pdf}
	\caption{Surface de cerne observée à 1,30 m en fonction de la surface de cerne prédite par le modèle.}
	\label{surf130_obs_vs_surf130_pred_sans_derive}
\end{figure}

\begin{figure}[!h]
	\centering	
	\includegraphics[width=\textwidth]{volume_obs_vs_volume_pred_sans_derive.pdf}
	\caption{Volume du cerne observé en fonction du volume du cerne prédit par le modèle.}
	\label{volume_obs_vs_volume_pred_sans_derive}
\end{figure}

\FloatBarrier

La Fig. \ref{modelisation_croissance_arbre_moyen_1_sans_derive} présente l'évolution des résidus pour l'estimation du $D_{130}$ au cours du temps lorsqu'on modélise l'arbre moyen. Le calcul est fait en cumulant au cours du temps les diamètres à 1,30 m modélisés. Ici le modèle prend en entrée le vrai diamètre à 1,30 m de l'année précédente ($D130_{n-1}$) et la vraie largeur moyenne à 1,30 m des 5 derniers cernes de l'année précédente ($LC5_{n-1}$). Au travers de cette dernière variable le modèle prend partiellement en compte du climat et l'effet des éclaircies. Les résultats pour Ecouves montrent une tendance à l'augmentation de la surestimation du modèle avec l'âge. 

La Fig. \ref{modelisation_croissance_arbre_moyen_2_sans_derive} présente les résidus que l'on obtiendrait si on utilisait en entrée du modèle les sorties calculées à partir des années précédentes. Les résidus restent très importants et beaucoup d'arbres présentent toujours une forte dérive, mais le point positif est que le modèle ne présente plus de biais systématique (sur-estimation du $D_{130}$) à Ecouves.

\begin{figure}[!h]
	\centering	
	\includegraphics[width=\textwidth]{modelisation_croissance_arbre_moyen_1_sans_derive.pdf}
	\caption{Résidus du $D_{130}$ de l'arbre moyen modélisé à partir des variables mesurées.}
	\label{modelisation_croissance_arbre_moyen_1_sans_derive}
\end{figure}

\begin{figure}[!h]
	\centering	
	\includegraphics[width=\textwidth]{modelisation_croissance_arbre_moyen_2_sans_derive.pdf}
	\caption{Résidus de la croissance en $D_{130}$ de l'arbre moyen modélisé à partir de variables prédites (processus itératif).}
	\label{modelisation_croissance_arbre_moyen_2_sans_derive}
\end{figure}

\FloatBarrier


\subsection{Modélisation des incréments annuels normés de la surface du cerne}

\label{Modele_Deleuze_normé}

Le but de cette approche de modélisation sur les profils d'incréments annuels normés est une implémentation du modèle de répartition dans le simulateur de croissance \textit{SimCop}. En effet dans le simulateur la quantité de matière à allouer chaque année est donnée ($BI$ pour \textit{bole increment}) et dépend de la croissance du houppier. Nous ne devons pas recalculer cette quantité mais uniquement la distribuer le long de la tige. C'est pourquoi nous avons décidé dans cette partie de travailler uniquement sur la forme du profil d'incrément de façon à ce que le volume du cerne soit égal à 1.

L'incrément annuel normé à une hauteur donnée est le rapport de l'incrément en surface par l'intégrale des incréments sur toute la hauteur de l'arbre, autrement dit par le volume du cerne divisé par la hauteur totale. Le volume du cerne a été calculé par intégration du profil sur toute la hauteur de l'arbre en interpolant linéairement les incréments en surface mesurés. L'incrément en surface a été fixé à 0 au niveau de l'apex et extrapolé à la hauteur 0 à partir des incréments mesurés sur les deux premières hauteurs de mesure (à 30 et 80 cm de hauteur pour Ecouves, à 1,30 m et environ 3 m pour Amance).

\subsubsection{Ajustement d'un modèle dérivé de \cite{Deleuze2002}}

Une version modifiée du modèle de \cite{Deleuze2002} a été utilisée pour prédire le profil d'incrément normé $\Gamma$ en fonction de la distance relative à l'apex $x_r$ (équation \ref{ShX_norme}).
\begin{equation}
\Gamma(x_r) = \frac{ \Delta G(x_r)}{I_G} = x_r \frac{P_r \sinh(z  x_r) + x_r  P_f \sinh(z(1-x_r))}{\sinh(z)}
\label{ShX_norme}
\end{equation}
avec $\Delta G(x_r)$ l'incrément ligneux annuel en surface à la distance relative $x_r$ à l'apex , $1-x_r$ étant égal à la hauteur relative. $I_G$ est l'intégrale de $\Delta G$ sur toute la hauteur de l'arbre.

\subsubsection{Modélisation des profils d'incréments annuels à partir de variables "arbre"}

La Fig. \ref{pairs_Pr_Pf_z_vs_var_arbre_norme} montre les corrélations existant entre les 3 paramètres du modèle et des variables intrinsèques à l'arbre.  Dans \textit{SimCop}, l'âge, la hauteur totale et le volume du cerne sont connus à l'appel de la fonction de répartition de l'accroissement courant. 
Il serait possible d'utiliser également le $D_{130}$ de l'année $n - 1$ mais, pour éviter tout risque de dérive par  accumulation des erreurs, nous avons préféré reprendre les variables  $D\_volume$, $H/D\_volume$ et $LC5\_volume$ (définies page \pageref{definitions_variables_volume}) dont le calcul est complètement indépendant de la forme du profil de cerne.
Pour ajuster les modèles, nous n'avons conservé ici que les profils des cernes d'âge supérieur à 10 ans. Un cerne très atypique (cerne 2010 de l'arbre 10 d'Ecouves) a également été éliminé.

\begin{figure}[!h]
	\centering
	\includegraphics[width=\textwidth]{pairs_Pr_Pf_z_vs_var_arbre_vol_norme}
	\caption{Corrélations entre différentes variables niveau "arbre" et les paramètres $P_{r}$, $P_{f}$ et $z$ du modèle \ref{ShX_norme} (sites d'Ecouves et d'Amance, cernes d'âge supérieur à 10ans).}
	\label{pairs_Pr_Pf_z_vs_var_arbre_norme}
\end{figure}

Nous avons choisi de modéliser les trois paramètres en fonction de l'âge.
Une première analyse des résidus à ce stade a montré l'existence d'un biais avec le rapport $H/D$ et la largeur des cinq derniers cernes pour les paramètres $P_{f}$ et $z$.
Bien que des relations assez fortes existent entre ces paramètres et la variable $LC5\_volume$, les modèles finalement adoptés utilisent plutôt le rapport $H/D\_volume$. Ils sont décrits par les équations \ref{equa_modelisation_z_norme} à \ref{equa_modelisation_Pr_norme}.

\begin{equation} 
	z = a_z \times \left(  1-\exp \left[ -H/D\_vol_n \times \frac{Age_n - c_z}{b_z} \right]  \right)
	\label{equa_modelisation_z_norme}
\end{equation}

avec $a_z = 5,524$, $b_z = 1812$, $c_z = 3,065$.

\begin{equation} 
	P_f = a_f \times H/D\_vol_n \times \left( 1 - \exp \left[ -\frac{Age_n}{b_f}  \right] \right)
	\label{equa_modelisation_Pf_norme}
\end{equation} 

avec $a_f = 0,4327$, $b_f = 17,63$.

\begin{equation} 
	% "P_r = ar * (age/br + 1/(1 - exp(-age/br)))"
	P_r = a_r \times \frac{Age_n}{b_r}  + \frac{1}{1-\exp  \left[-\frac{Age_n}{b_r} \right]}
	\label{equa_modelisation_Pr_norme}
\end{equation} 

avec $a_r = 0,6404$, $b_r = 23,28$.

$Age_n$ est l'âge (en année) et$H/D\_vol_n$ le rapport hauteur sur diamètre (sans dimension) défini page \pageref{definitions_variables_volume} pour l'année courante.

Les figures \ref{modelisation_z_norme} à \ref{modelisation_Pr_residus_norme} montrent graphiquement les ajustements obtenus et l'analyse des résidus correspondants. Le comportement est globalement satisfaisant même si quelques imperfections apparaissent pour les plus jeunes arbres et les faibles vitesses de croissance.

\begin{figure}[!h]
	\centering	
	\includegraphics[width=9cm, page=3]{modelisation_Pr_Pf_z_norme}
	\caption{Modélisation du paramètre $z$ en fonctionde l'âge de l'année $n-1$. Les trois courbes bleues sont obtenues avec des valeurs différentes du rapport $H/D$ (moyenne, premier et dernier vingtile).}
	\label{modelisation_z_norme}
\end{figure}

\begin{figure}[!h]
	\centering
	\includegraphics[width=\textwidth, page=3]{modelisation_Pr_Pf_z_residus_norme}
	\caption{Résidus du modèle en fonction de variables "arbre".}
	\label{modelisation_z_residus_norme}
\end{figure}

\begin{figure}[!h]
	\centering
	\includegraphics[width=9cm, page=2]{modelisation_Pr_Pf_z_norme}
	\caption{Modélisation du paramètre $P_f$ en fonction de l'âge de l'année $n-1$. Les trois courbes bleues sont obtenues avec des valeurs différentes du rapport $H/D$ (moyenne, premier et dernier vingtile).}
	\label{modelisation_Pf_norme}
\end{figure}

\begin{figure}[!h]
	\centering
	\includegraphics[width=\textwidth, page=2]{modelisation_Pr_Pf_z_residus_norme}
	\caption{Résidus du modèle en fonction de variables "arbre".}
	\label{modelisation_Pf_residus_norme}
\end{figure}

\begin{figure}[!h]
	\centering
	\includegraphics[width=9cm, page=1]{modelisation_Pr_Pf_z_norme}
	\caption{Modélisation du paramètre $P_r$ en fonction de l'âge de l'année $n-1$.}
	\label{modelisation_Pr_norme}
\end{figure}

\begin{figure}[!h]
	\centering
	\includegraphics[width=\textwidth, page=1]{modelisation_Pr_Pf_z_residus_norme}
	\caption{Résidus du modèle en fonction de variables "arbre".}
	\label{modelisation_Pr_residus_norme}
\end{figure}

\FloatBarrier
\bigskip

Nous avons ensuite ajusté le modèle global (équations \ref{equa_modele_global_norme}) en recherchant des effet site, traitement et arbre sur les paramètres du modèle avec la fonction \textit{nlsList}. La figure \ref{nlsList_effet_site_norme} montre qu'un effet site existe sur tous les paramètres du modèle.  Indépendamment de l'effet site, un effet traitement a été observé pour $a_r$ et $b_r$ sur les deux sites (Fig. \ref{nlsList_effet_traitement_norme}) et un effet arbre dans traitement sur $a_r$, $b_r$ et $a_z$.

Nous avons choisi de rentrer ces effets arbre, traitement et site en effets aléatoires dans le modèle final.

\begin{equation} 
	\begin{array}{l}
	z = a_z \times \left(  1-\exp \left[ -H/D\_vol_n \times \frac{Age_n - c_z}{b_z} \right]  \right)\\
	P_f = a_f \times H/D\_vol_n \times \left( 1 - \exp \left[ -\frac{Age_n}{b_f}  \right] \right)\\
	P_r = a_r \times \frac{Age_n}{b_r}  + \frac{1}{1-\exp  \left[-\frac{Age_n}{b_r} \right]}\\
	\\
	\Gamma(x_r) = x_r \frac{P_r \sinh(z  x_r) + x_r  P_f \sinh(z(1-x_r))}{\sinh(z)}\\
	\end{array}
	\label{equa_modele_global_norme}
\end{equation} 

où $\Gamma(x_r)$ est l'incrément en surface normé à la distance relative $x_r$ à l'apex,
$Age_n$ est l'âge (en année) et $H/D\_vol_n$ le rapport hauteur sur diamètre (sans dimension) défini page \pageref{definitions_variables_volume} pour l'année courante.\\

Les valeurs estimées des paramètres (effets fixes) du modèle sont les suivantes :
$a_r = 0,61528$, $b_r = 25,241$, $a_f = 0,42728$, $b_f = 21,151$, $a_z = 6,0926$, $b_z = 5373,2$, $c_z = -24,448$.

\begin{figure}[!h]
	\centering	
	\includegraphics[width=\textwidth, page=1]{nlsList_norme}
	\caption{Recherche d'un effet site sur les paramètres du modèle.}
	\label{nlsList_effet_site_norme}
\end{figure}

\begin{figure}[!h]
	\centering	
	\includegraphics[width=\textwidth, page=2]{nlsList_norme}
	\caption{Recherche d'un effet traitement sur les paramètres du modèle.}
	\label{nlsList_effet_traitement_norme}
\end{figure}

Après la procédure d'ajustement du modèle les résultats ne sont plus forcément normés. Pour obtenir l'incrément en surface normé, il est donc nécessaire de diviser les valeurs estimées à chaque hauteur par l'intégrale des valeurs estimées sur tout le profil.

\begin{equation} 
	\begin{array}{l}
	\Gamma_r(x_r) = \frac{\Gamma(x_r)}{I_{\Gamma}}\\
	\end{array}
	\label{equa_renormage}
\end{equation}

où $\Gamma_r(x_r)$ est l'incrément en surface normé à la distance relative $x_r$ à l'apex, $\Gamma(x_r)$ est la valeur estimée par le modèle \ref{equa_modele_global_norme} et $I_{\Gamma}$ est le résultat de l'intégration numérique de $\Gamma$ sur l'intervalle $[0, 1]$.

L'estimation de l'incrément absolu en surface peut finalement être obtenue par l'équation \ref{equa_denormage} :

\begin{equation} 
	\begin{array}{l}
	\Delta G(x_r) = \Gamma_r(x_r) \times \frac{V_{cerne}}{H}\\
	\end{array}
	\label{equa_denormage}
\end{equation}

où $\Delta G(x_r)$ est l'incrément en surface à la distance relative $x_r$ à l'apex, $\Gamma_r(x_r)$ est  l'incrément en surface normé, $V_{cerne}$ est le volume du cerne et $H$ la hauteur totale.

Les résidus du modèle complet (Fig. \ref{residus_absolus_modele_norme}) ne présentent pas de tendance marquée avec les autres variables descriptives. On note néanmoins que le modèle tend à sous-estimer les surfaces des cernes les plus larges de l'empattement (mesures d'Ecouves à 30 cm de hauteur). Ce défaut s'aggrave sensiblement lorsque l'on ne considère que la partie fixe du modèle (Fig. \ref{residus_absolus_modele_partie_fixe_norme}).

\begin{figure}[!h]
	\centering	
	\includegraphics[width=\textwidth]{residus_absolus_modele_mixte_arbre_norme}
	\caption{Résidus du modèle global (équations \ref{equa_modele_global_norme} et \ref{equa_denormage}) avec effets aléatoires site, traitement et arbre sur les paramètres du modèle.}
	\label{residus_absolus_modele_norme}
\end{figure}
\begin{figure}[!h]
	\centering	
	\includegraphics[width=\textwidth]{residus_absolus_modele_mixte_arbre_partie_fixe_norme}
	\caption{Résidus de la partie fixe du modèle global (équations \ref{equa_modele_global_norme} et \ref{equa_denormage}).}
	\label{residus_absolus_modele_partie_fixe_norme}
\end{figure}

% Question à poser à Fleur :
% Pourquoi les résidus sont plus forts avec le modèle complet pour certains arbres ??
% 	Cf surf130_obs_vs_surf130_pred_modele_mixte_arbre_norme.pdf
% 	et surf130_obs_vs_surf130_pred_modele_mixte_arbre_partie_fixe_norme.pdf

Au delà de l'empattement, les résultats sont satisfaisants. 
La Fig. \ref{surf130_obs_vs_surf130_pred_norme} représente les valeurs observées en fonction des valeurs prédites pour la surface de cerne à 1,30 m.
La RMSE pour la surface de cerne à 1,30 m est de $4,05 \times 10^{-4}$ m$^{2}$ et la RMSE relative correspondante est de 24 \%.
La RMSE à Ecouves est de $4,41 \times 10^{-4}$ m$^{2}$ tandis qu'à Amance elle est de $2,82 \times 10^{-4}$ m$^{2}$

L'erreur moyenne est de $0,17 \times 10^{-4}$ m$^{2}$ à Ecouves (1 \%) tandis qu'elle n'est que de $0,02 \times 10^{-4}$ m$^{2}$ à Amance (0,1 \%). Le modèle sous-estime donc légèrement la surface de cerne à 1,30 m.

\begin{figure}[!h]
	\centering	
	\includegraphics[width=\textwidth]{surf130_obs_vs_surf130_pred_norme.pdf}
	\caption{Surface de cerne observée à 1,30 m en fonction de la surface de cerne prédite par le modèle.}
	\label{surf130_obs_vs_surf130_pred_norme}
\end{figure}

La Fig. \ref{modelisation_croissance_arbre_moyen_norme} présente l'évolution des résidus pour l'estimation du $D_{130}$ au cours du temps lorsqu'on modélise l'arbre moyen (en utilisant uniquement la partie fixe du modèle). Le calcul est fait en cumulant au cours du temps les largeurs de cernes prédites à 1,30 m. Ici le modèle prend en entrée la hauteur totale et le volume du cerne de l'année courante conformément aux conditions d'utilisation dans \textit{SimCop}. Le volume du cerne reflète évidemment les conditions de croissance, y compris le climat, ce qui explique que les résidus soient ici sensiblement inférieurs à ceux des modèles précédents (Fig. \ref{modelisation_croissance_arbre_moyen_1_sans_derive} et Fig. \ref{modelisation_croissance_arbre_moyen_2_sans_derive}). La figure \ref{diametre_mesures-estimes_modele_mixte_arbre_partie_fixe_norme} permet de vérifier l'absence de biais notable dans les $D_{130}$ prédits.

\begin{figure}[!h]
	\centering	
	\includegraphics[width=\textwidth]{modelisation_croissance_arbre_moyen_norme}
	\caption{Résidus du $D_{130}$ de l'arbre moyen modélisé à partir des variables mesurées.}
	\label{modelisation_croissance_arbre_moyen_norme}
\end{figure}

\begin{figure}[!h]
	\centering	
	\includegraphics[width=9cm]{diametre_mesures-estimes_modele_mixte_arbre_partie_fixe_norme}
	\caption{Diamètres à 1,30 m mesurés et prédits par le modèle à l'âge final (39 ans pour Amance, 63 ans pour Ecouves).}
	\label{diametre_mesures-estimes_modele_mixte_arbre_partie_fixe_norme}
\end{figure}

L'annexe \ref{Annexe_profils_normes} montre l'évolution avec l'âge des profils obtenus avec la partie fixe du modèle pour les arbres des deux dispositifs. Au stade juvénile, les profils obtenus sont quasiment coniques et évoluent progressivement vers une courbe en \textit{S}. On notera que l'effet du statut social est de plus en plus perceptible avec le temps. Dans les dernières années de croissance, les arbres dominés présentent une partie médiane nettement oblique là où le profil des arbres dominants est presque vertical.

\FloatBarrier

Pour l'année de sécheresse de 1976 nous avons tracé les résidus en fonction de la hauteur relative dans l'arbre (Fig. \ref{modelisation_globale_residus_vs_hauteur_relative_1976_modele_mixte_arbre_partie_fixe_norme}) pour vérifier la concordance avec des observations de la littérature qui relatent qu'en cas de stress intense les profils d'incréments annuels sont plus réduit dans le bas des tiges (avec même parfois pour certaines essences comme le Hêtre la disparition complète du cerne) que dans le haut (allocation prioritaire en haut).
La Fig. \ref{residus-abs_vs_annee_par_hauteur_modele_mixte_arbre_partie_fixe_norme} montre cependant que les résidus du modèle dans le bas de la tige (1,30 m) ne sont pas forcément plus faibles les années supposées de sécheresse.

\begin{figure}[!h]
	\centering	
	\includegraphics[width=\textwidth]{modelisation_globale_residus_vs_hauteur_relative_1976_modele_mixte_arbre_partie_fixe_norme}
	\caption{Résidus en fonction de la hauteur relative pour l'année de sécheresse de 1976.}
	\label{modelisation_globale_residus_vs_hauteur_relative_1976_modele_mixte_arbre_partie_fixe_norme}
\end{figure}

\begin{figure}[!h]
	\centering	
	\includegraphics[width=\textwidth]{residus-abs_vs_annee_par_hauteur_modele_mixte_arbre_partie_fixe_norme}
	\caption{Résidus en surface de cerne absolue (moyenne et intervalle de confiance) du modèle global (Équation \ref{equa_denormage}, partie fixe du modèle mixte) en fonction de l'année calendaire. Les lignes verticales grises en pointillés correspondent aux dates des éclaircies. Les lignes verticales rouges en pointillés correspondent à des années de sécheresse (à vérifier). Les données sont présentées par traitement sylvicole et pour trois hauteurs : 30 cm, 1,30 m et la rondelle située juste sous la base du houppier estimée. Seul le site d'Ecouves est représenté car sur le site d'Amance il n'y a pas de rondelles prélevées sous 1,30 m.}
	\label{residus-abs_vs_annee_par_hauteur_modele_mixte_arbre_partie_fixe_norme}
\end{figure}

\begin{figure}[!h]
	\centering	
	\includegraphics[width=\textwidth]{residus-normes_vs_annee_par_hauteur_modele_mixte_arbre_partie_fixe_norme}
	\caption{Résidus en surface de cerne normée (moyenne et intervalle de confiance) du modèle global (Équation \ref{equa_renormage}, partie fixe du modèle mixte) en fonction de l'année calendaire. Les lignes verticales grises en pointillés correspondent aux dates des éclaircies. Les lignes verticales rouges en pointillés correspondent à des années de sécheresse (à vérifier). Les données sont présentées par traitement sylvicole et pour trois hauteurs : 30 cm, 1,30 m et la rondelle située juste sous la base du houppier estimée. Seul le site d'Ecouves est représenté car sur le site d'Amance il n'y a pas de rondelles prélevées sous 1,30 m.}
	\label{residus-normes_vs_annee_par_hauteur_modele_mixte_arbre_partie_fixe_norme}
\end{figure}

\FloatBarrier
\subsubsection{Implémentation dans \textit{SimCop}}

Le modèle a été implémenté dans une version de travail du simulateur \textit{SimCop} sous \textit{Capsis} en laissant la possibilité à l'utilisateur de conserver la méthode d'origine (allocation basée sur la loi de Pressler) ou d'utiliser le nouveau modèle (basé sur les équations \ref{equa_modele_global_norme} et \ref{equa_denormage}).

Une première simulation d'une plantation à 3 m $\times$ 3 m laissée en croissance libre pendant 60 ans a été réalisée avec les deux options. En ce qui concerne le diamètre à 1,30 m, les résultats obtenus (Fig. \ref{simcop_capsis}) montrent que le diamètre quadratique moyen est identique dans les deux cas mais que le diamètre des arbres dominants est légèrement plus faible avec le nouveau modèle. La figure \ref{simcop_profil} montre que les profils de tige obtenus sont sensiblement différents, notamment au niveau de l'empattement.

\begin{figure}[!h]
	\centering
	\includegraphics[width=\textwidth]{simcop_capsis}
	\caption{Simulation sous capsis du même scénario jusqu'à 60 ans avec la méthode d'origine (violet) et le nouveau modèle (rouge) d'allocation de l'incrément annuel.}
	\label{simcop_capsis}
\end{figure}

\begin{figure}[!h]
	\begin{tabular}{cc} 
		\includegraphics[width=.5\textwidth]{simcop_profil_ori}
		\includegraphics[width=.5\textwidth]{simcop_profil_nou}\\
	\end{tabular}
	\caption{Profil de tige simulé à 60 ans avec la méthode d'origine (à gauche) et le nouveau modèle (à droite) sur le même arbre pour lequel le diamètre à 1,30 m est presque identique (23,8 cm à gauche, 23,7 cm à droite).}
	\label{simcop_profil}
\end{figure}

Pour le même scénario nous avons comparé les hauteurs de base de houppier prédites à 60 ans avec le profil de cerne en utilisant l'équation \ref{pred_base_houppier} (page \pageref{pred_base_houppier}) avec les hauteurs de base de houppier calculées par le moteur de Simcop.
La figure \ref{simcop_hbh_est} montre qu'il existe une certaine cohérence entre les deux variables mais que la hauteur obtenue avec le profil de cerne est sensiblement plus basse.


Un travail approfondi de simulation sera nécessaire pour tester le modèle et analyser son comportement dans différentes situations de croissance avant de décider de valider ou non cette modification.

\begin{figure}[!h]
	\centering
	\includegraphics[width=.75\textwidth]{simcop_hbh_est}
	\caption{Relation entre hauteur de base de houppier prédite par l'équation \ref{pred_base_houppier} (HBH\_est) et hauteur de base de houppier calculée par Simcop (HBH) à 60 ans.}
	\label{simcop_hbh_est}
\end{figure}

\FloatBarrier
\section{Perspectives}

Nous souhaitons - entre autres- utiliser ce jeu de données pour tester si les formes de tiges générées par le simulateur \textit{SimCop} sont réalistes par rapport à la sylviculture pratiquée. L'historique du peuplement est connue (dates des interventions et nombre d'arbres supprimés à chaque date). En revanche les éclaircies ne sont pas spatialisées. L'idée serait de simuler plusieurs scénarios en accord avec l'intensité des éclaircies observées pour étudier l'étendue des variations générées et comparer les résultats obtenus aux formes de tiges réelles. Ce travail serait fait de façon coordonnée avec un autre travail en cours sur la validité des sorties des simulateurs de croissance existants du Douglas (post-doc en cours co-financé par la chaire ONF-APT et la DGER).\\

Par ailleurs, nous avons travaillé ici sur un jeu de données certes très détaillées (mesures de largeurs de cernes à plusieurs hauteurs) mais sur un nombre d'arbres très faible. Par la suite, il serait possible de travailler sur un grand nombre d'arbres (par exemple base de données \textit{EMERGE} mais pas seulement) de façon à couvrir une grande diversité de situations. En revanche sur ces arbres seul un profil de tige à une date donnée sera disponible. La question est de savoir si pour ces arbres on aura assez d'information disponible sur la gestion passée du peuplement pour pouvoir en tirer des conclusions.\\

L'étape ultérieure de ce travail, à court terme (démarrage prévu en 2014), sera de travailler sur les effets des traitements sylvicoles sur les incréments annuels de croissance le long de la tige de façon à pouvoir \textit{in fine} développer un modèle dynamique de croissance étant à même de bien décrire les changements de défilement et de forme des tiges liés à la sylviculture. 


\chapter{Comprendre et modéliser la distribution de l'aubier des tiges de Douglas en France}

Ce travail a été initié lors du stage de M2 de Marwen Fatnassi (Master 2 "Ingénierie Mathématique et Outils Informatiques") financé par la convention Modelfor 2012-2015. Seul le jeu d'Ecouves contient des mesures de largeurs d'aubier à plusieurs hauteurs. Le jeu d'Amance ne dispose pas de telles mesures.

\section{Revue de la littérature}

\subsection{Les différents modèles de distribution de l'aubier et du duramen}

\subsection{Les effets des traitements sylvicoles sur la distribution de l'aubier et du duramen}

\section{Résultats}

\subsection{Analyse préliminaire}

Nous avons choisi de regrouper les deux traitements éclaircis "B+C" car ils n'étaient pas suffisamment contrastés.
 
Le tableau \ref{tableau_caracteristiques} montre que l'effet du statut social (et donc du diamètre) est prépondérant par rapport à l'effet du traitement sylvicole. Entre un dominé et un dominant des traitements éclaircis le rayon du duramen augmente de 6,8 cm tandis que la largeur d'aubier n'augmente que de 2,6 cm. Un arbre dominant a donc à la fois plus de duramen et plus d'aubier. Les arbres des peuplements éclaircis ont également plus de duramen et plus d'aubier que ceux du peuplement témoin ("A") mais dans une moindre mesure. Ce tableau reste difficile à interpréter dans la mesure où les dominants des peuplements éclaircis ont des diamètres supérieurs à ceux du traitement témoin. Des analyses statistiques plus poussées via l'utilisation de modèles multivariés devraient nous permettre d'y voir plus clair. 

\begin{table}[h]
	\caption{Moyennes à 6 m du rayon du duramen, de la largeur d'aubier et de la proportion d'aubier en surface en fonction du traitement sylvicole ("A" = témoin vs. "B+C") et du statut social.}
	\centering 
		\begin{tabular}{lllll} 
			\hline
			Traitement & Statut social & Rayon du & Largeur & Proportion d'aubier \\
			& & duramen (cm) & d'aubier (cm) & en surface  (\%) \\
			\hline 
			A & Dominant & 15,1 & 3,6 & 34,8 \\
			B+C & Dominant & 16,1 & 5,0 & 41,6 \\			
			\hline 
			A & Co-dominant & 11,3 & 2,7 & 35,4 \\
			B+C & Co-dominant & 13,0 & 3,7 & 39,5 \\			
			\hline
			A & Dominé & 9,4 & 2,1 & 33,2 \\
			B+C & Dominé & 9,3 & 2,4 & 36,3 \\			
			\hline
		\end{tabular}
	\label{tableau_caracteristiques}
\end{table}

Le but de cette étude est de développer un modèle de défilement de l'aubier et du duramen. Les différentes variables sont étudiées de façon à choisir celles qui seront les plus pertinentes à modéliser.

\bigskip 

Liste des variables potentielles à modéliser :
\begin{itemize}
\item Rayon du duramen
\item Largeur d'aubier
\item Proportion d'aubier en largeur (ou de duramen au choix)
\item Surface du duramen
\item Surface d'aubier
\item Proportion d'aubier en surface (ou de duramen au choix)
\item Nombre de cernes du duramen
\item Nombre de cernes d'aubier
\item Proportion d'aubier en nombre de cernes (ou de duramen au choix)
\end{itemize}

\bigskip 

Les profils par arbre des différentes variables sont représentés en annexe \ref{Appendix_LA}. 

Les Fig. \ref{variables_vs_hauteur_par_statut} et \ref{variables_vs_hauteur_par_traitement} représentent les variations des mêmes variables en fonction de la hauteur dans l'arbre (entre 1,30 m et la base du houppier uniquement ici) par statut social et par traitement sylvicole, respectivement. Nous remarquons notamment qu'il semble y avoir un fort effet du statut social sur ces variables alors que l'effet du traitement sylvicole semble moins évident à première vue.

\bigskip 

Liste des variables explicatives au niveau "arbre" :
\begin{itemize}
\item Hauteur totale de l'arbre
\item Hauteur de la base du houppier
\item Hauteur de la première branche verte
\item $C_{130}$
\item Statut social
\item Surface projetée du houppier
\item Largeur moyenne de cerne à 1,30 m
\item Largeur moyenne à 1,30 m des 5 derniers cernes (2006 à 2010)
\item $H/D$
\item Longueur relative du houppier $LHP/HT$
\end{itemize}

\bigskip 

Liste des variables explicatives au niveau "peuplement" :
\begin{itemize}
\item Traitement sylvicole
\end{itemize}

\bigskip 

A titre de comparaison avec la littérature, nous trouvons une diminution de la surface d'aubier entre 1,30 m et la base du houppier de 61\% en moyenne. \cite{Waring1982} sur 16 Douglas de 24 m de haut en moyenne et avec des couronnes vivantes représentant 48\% de la hauteur de l'arbre en moyenne avaient trouvé une décroissance de 36\% seulement.

\begin{figure}[!h]
	\centering
	\includegraphics[width=\textwidth]{variables_vs_hauteur_par_statut.pdf}	
	\caption{Représentation des différentes variables en fonction de la hauteur dans l'arbre par statut social. La ligne verticale grise en pointillés indique la hauteur de 1,30 m.}
	\label{variables_vs_hauteur_par_statut}
\end{figure}

\begin{figure}[!h]
	\centering
	\includegraphics[width=\textwidth]{variables_vs_hauteur_par_traitement.pdf}	
	\caption{Représentation des différentes variables en fonction de la hauteur dans l'arbre par traitement sylvicole. La ligne verticale grise en pointillés indique la hauteur de 1,30 m.}
	\label{variables_vs_hauteur_par_traitement}
\end{figure}

\FloatBarrier


\subsection{Modélisation de la largeur d'aubier (modèles 1 et 4)}

\subsubsection{Modèle 1}

Le modèle 1 est basé sur l'observation que la largeur d'aubier est relativement constante entre 1,30 m et la base du houppier. Par conséquent une largeur d'aubier constante est modélisée dans cette partie de l'arbre. A l'empattement c'est un segment de droite qui est ajusté à partir des 3 rondelles disponibles dans cette partie (0,30 m, 0,80 m et 1,30 m).  

\paragraph{A l'empattement (entre le sol et 1,30 m) :}

L'annexe \ref{largeur_aubier_empattement} montre l'évolution de la largeur d'aubier à l'empattement. Les hauteurs de rondelles sont fixes dans cette partie de l'arbre : 0,30 m, 0,80 m et 1,30 m. Comme première approche\footnote{Il sera toujours possible de raffiner après en modélisant une courbure.} nous avons choisi d'ajuster une droite à partir des 3 points et de modéliser la pente $\alpha$ de cette droite à partir des variables "arbre" et "peuplement" disponibles. Nous avons uniquement besoin de modéliser la pente car à 1,30 m on fixe la largeur d'aubier à la valeur constante que nous aurons modélisée dans le tronc entre 1,30 m et la base du houppier. 
Un graphique des corrélations montre que les variables les plus corrélées à la pente $\alpha$ sont dans l'ordre : la largeur moyenne des 5 derniers cernes ($r=-0,50$, $\text{p-value}=0,00760$), la surface de la projection du houppier ($r=-0,46$, $\text{p-value}=0,01619$), la longueur relative du houppier ($r=-0,43$, $\text{p-value}=0,02347$), puis le $C_{130}$ ($r=-0,42$, $\text{p-value}=0,02913$). Par souci d'obtenir un modèle facilement applicable nous avons retenu un modèle incluant le $C_{130}$ (Équation \ref{pente_empattement}). Il y a de plus un effet du traitement sylvicole à la fois sur les valeurs de $a$ et de $b$ ci-dessous (Fig. \ref{LA_vs_C130_empattement_par_traitement} et Tableau \ref{table_estimation_parametres_empattement}). 
Assez curieusement il apparaît que le traitement B se distingue du A mais pas le C qui en théorie devrait être plus éloigné encore. Pour l'application des modèles nous avons décidé de faire un groupe contenant le traitements "éclaircis" B et C sans les distinguer. Sur la Fig. \ref{LA_vs_C130_empattement_par_traitement}, nous remarquons qu'il est possible d'avoir des pentes positives (c'est-à-dire une diminution de la largeur d'aubier à l'empattement) pour les petits $C_{130}$ dans les peuplements éclaircis et à l'inverse pour les grands $C_{130}$ dans le peuplement témoin.

\begin{equation}
	\boxed{ 
	\alpha = a_{traitement} + b_{traitement} \times C_{130}(cm) 
	\label{pente_empattement}
	}
\end{equation}

\begin{table}[h]
	\caption{Estimation des paramètres du modèle correspondant à l'équation \ref{pente_empattement}.}
	\centering 
		\begin{tabular}{rrrr} 
			\hline
			Traitement & $n$ & $a$ & $b$ \\
			\hline 
			Tous & 27 & 0,01064 & -0,00014 \\
			\hline 
			A & 9 & -0,01643 & 0,00013 \\
			B & 9 & 0,02676 & -0,00031 \\ 
			C & 9 & -0,00098 & -4e-05 \\
			\hline
			B+C & 18 & 0,01574 & -0,00019 \\
			\hline
		\end{tabular}
	\label{table_estimation_parametres_empattement}
\end{table}

\begin{figure}[!h]
	\centering
	\includegraphics[width=9cm]{LA_vs_C130_empattement_par_traitement.pdf}	
	\caption{Relation entre le diamètre à 1,30 m et la pente de la largeur d'aubier à l'empattement (entre le sol et 1,30 m) représentée par traitement sylvicole (Équation \ref{pente_empattement} avec les valeurs de paramètres du Tableau \ref{table_estimation_parametres_empattement}).}
	\label{LA_vs_C130_empattement_par_traitement}
\end{figure}

\FloatBarrier

\paragraph{Entre 1,30 m et la base du houppier :}

Seules les rondelles situées entre 1,30 m et la base du houppier de chaque arbre ont été conservées pour les analyses. Un arbre inférieur à 1,30 m n'aura de toute façon pas de valeurs de largeur d'aubier ($LA$) dans cette zone. Pour un arbre de 1,30 m exactement on aura $C_{130} = 0$ et $LA = 0$ à 1,30 m et pour un arbre $> 1,30$ m on aura $LA > 0$ partout au-dessus de 1,30 m car il y a toujours de l'aubier (ce n'est pas le cas du duramen). En fait, plus exactement le modèle ne pourra être appliqué raisonnablement que si la hauteur de la base du houppier est au-dessus de 1,30 m car il est calibré dans cette zone (au-dessus de l'empattement et en-dessous de la base du houppier) sur des arbres matures de 63 ans.

La largeur d'aubier qui paraissait assez constante le long du tronc sur les graphiques au niveau arbre (Annexe \ref{Appendix_LA}) présente en fait un schéma de forme en cuvette qui se répète pour chaque arbre. Les variations observées restent cependant relativement faibles. Des régressions linéaires par arbre en fonction de la hauteur dans l'arbre montrent que l'effet de la hauteur est souvent non significatif et que les valeurs de pente sont très proches de 0, comprises entre $-6,8\times10^{-4}$ et $3,3\times10^{-4}$. Comme déjà montré dans d'autres études \citep[e.g.,][]{Longuetaud2006}, le \textit{pipe model} de \cite{shinozaki1964a,shinozaki1964b} n'est pas valide sous le houppier car s'il l'était alors la surface d'aubier devrait être constante ce qui n'est pas le cas.

Dans un premier temps, et pour simplifier le problème, nous considérons une valeur de largeur d'aubier constante, notée $LA$, dans cette partie de l'arbre et nous cherchons à la modéliser. La variable à modéliser est la largeur moyenne d'aubier observée entre 1,30 m et la base du houppier. La Figure \ref{pairs_LA_moy_vs_var_arbre} illustre les corrélations observées entre les variables mesurées au niveau de l'arbre et la largeur d'aubier moyenne. La meilleure variable pour expliquer la largeur d'aubier sous le houppier est la largeur de cerne moyenne à 1,30 m ($r=0,91$, $\text{p-value}=5,734e-11$), suivie de la circonférence à 1,30 m ($r=0,89$, $\text{p-value}=3,257e-10$). Nous retiendrons plutôt cette dernière variable car elle est beaucoup plus facile à obtenir sur le terrain. Un premier ajustement a montré que le paramètre d'ordonnée à l'origine était non significatif. 

\begin{figure}[!h]
	\centering
	\includegraphics[width=\textwidth]{pairs_LA_moy_vs_var_arbre.pdf}	
	\caption{Corrélations entre les différentes variables niveau "arbre" et la largeur d'aubier moyenne entre 1,30 m et la base du houppier.}
	\label{pairs_LA_moy_vs_var_arbre}
\end{figure}

Une analyse des résidus du modèle \ref{modele_2_LA_sous_houppier_sans_traitement} est présentée en Annexe \ref{largeur_aubier_tronc_residus}. Il s'avère qu'il reste un effet du traitement sylvicole visible notamment sur les résidus du modèle et statistiquement significatif quand introduit dans le modèle sur le paramètre de pente. En effet, la largeur d'aubier est surestimée dans le peuplement témoin A alors qu'elle est sous-estimée dans les peuplements éclaircis B et C. L'équation finalement retenue est donc l'équation \ref{modele_2_LA_sous_houppier}, sans ordonnée à l'origine mais avec un effet du traitement sylvicole sur la pente $a$. La pente obtenue est plus forte dans les deux peuplements éclaircis B et C en comparaison au témoin (A) (Tableau \ref{table_estimation_parametres} et Fig. \ref{LA_vs_C130_sous_houppier_par_traitement}). A âge et $C_{130}$ égaux, il y a plus d'aubier dans les arbres des peuplements éclaircis. La RMSE (en largeur moyenne de l'aubier) du modèle retenu dans cette partie de l'arbre (Équation \ref{modele_2_LA_sous_houppier}) est de 4,29 mm et la RMSE relative est de 12.6 \%. 

\begin{equation} 
	\boxed{
	LA (mm) = a \times C_{130}(cm) 
	\label{modele_2_LA_sous_houppier_sans_traitement}
	}
\end{equation}

\begin{equation} 
	\boxed{
	LA (mm) = a_{traitement} \times C_{130}(cm) 
	\label{modele_2_LA_sous_houppier}
	}
\end{equation}


\begin{table}[h]
	\caption{Estimation des paramètres du modèle correspondant à l'équation \ref{modele_2_LA_sous_houppier}.}
	\centering 
		\begin{tabular}{rrrr} 
			\hline
			Traitement & $n$ & $a$ \\
			\hline 
			Tous & 27 & 0,27054 \\
			\hline 
			A & 9 & 0,24110 \\
			B & 9 & 0,28823 \\ 
			C & 9 & 0,27571 \\
			\hline
			B+C & 18 & 0,28170 \\
			\hline
		\end{tabular}
	\label{table_estimation_parametres}
\end{table}

\begin{figure}[!h]
	\centering
	\includegraphics[width=9cm]{LA_vs_C130_sous_houppier_par_traitement.pdf}	
	\caption{Relation entre le diamètre à 1,30 m et la largeur moyenne d'aubier mesurée entre 1,30 m et la base du houppier représentée par traitement sylvicole (Équation \ref{modele_2_LA_sous_houppier} avec les valeurs de paramètres du Tableau \ref{table_estimation_parametres}).}
	\label{LA_vs_C130_sous_houppier_par_traitement}
\end{figure}


Le facteur traitement sylvicole pourrait peut-être être remplacé par d'autres variables intrinsèques à l'arbre. Un des meilleurs modèles obtenus à partir des variables "arbre" (\textit{stepwise} et suppression des variables \textit{ns}) est celui de l'équation \ref{modele_2_LA_sous_houppier_var_arbre} ($R^{2}=0,90$). La RMSE du modèle est de 3,32 mm et la RMSE relative est de 9,7 \%. Un seconde modèle possible est celui correspondant à l'équation \ref{modele_2_LA_sous_houppier_var_arbre_2}. La RMSE de ce modèle est de 3.86 mm et la RMSE relative est de 11,3\%. Les variables explicatives sont pour beaucoup fortement corrélées entre elles et plusieurs variantes de modèles sont possibles à définir en fonction des besoins d'applications, sans que les résultats ne changent très fortement. 

\begin{multline} 
	LA (mm) = a + b \times C_{130}(cm) + c \times Hauteur\_totale(m) + \\d \times Projection\_houppier(m^{2}) + e \times LC\_moyenne\_130(mm) 
	\label{modele_2_LA_sous_houppier_var_arbre}
\end{multline}

avec $a = 33,04961$, $b = -0,50382$, $c = 1,22046$, $d = 0,30573$ et $e = 24,45453$. 

\begin{multline} 
	LA (mm) = a \times Hauteur\_totale(m) + b \times Hauteur\_base\_houppier(m) + \\c \times LC\_moyenne\_130(mm) 
	\label{modele_2_LA_sous_houppier_var_arbre_2}
\end{multline}

avec $a=1,06369$, $b=-1,42968$ et $c=9,88221$.

\FloatBarrier

\paragraph{Dans le houppier}

L'annexe \ref{largeur_aubier_houppier} montre l'évolution de la largeur d'aubier dans le houppier. C'est dans cette partie de l'arbre que le \textit{pipe model} est censé être valide. Cependant pour le jeu de données d'Ecouves, nous ne disposons malheureusement pas du détails des branches, uniquement des projections de houppier. Dans la littérature la surface d'aubier à une hauteur donnée dans le houppier est souvent mise en relation avec la quantité de branches vivantes présentes au-dessus de ce niveau \citep[e.g.,][]{Schneider2011}, par exemple avec la somme des surfaces des sections de branches vivantes.

Il serait cependant possible de proposer une forme de modèle dans cette partie de l'arbre mais l'intérêt de le faire est limité.

\paragraph{Application du modèle en prédiction}

Nous avons fait le choix de développer un modèle facilement applicable prenant en entrée des variables facilement mesurables sur le terrain. Il est toujours possible de développer un modèle légèrement meilleur à partir de variables telles que la largeur moyenne de cernes à 1,30 m ou bien la largeur moyenne des 5 derniers cernes. Le modèle finalement développé prend donc en entrées la circonférence à 1,30 m et le traitement sylvicole "témoin" vs. "éclairci" et l'hypothèse est faite que la largeur d'aubier est constante entre 1,30 m et la base du houppier (ce qui n'est pas tout à fait vrai mais elle n'est pas très variable non plus). La Fig. \ref{Prediction_LA_appli_modele} montre les profils de largeur d'aubier obtenus sous le houppier pour des arbres de diamètre à 1,30 m 30 cm et 50 cm, respectivement, issus de peuplements éclaircis ou non. Nous retrouvons bien l'observation faite précédemment sur ce jeu de données, qu'à diamètre égal (même si on sait bien que les gros diamètres vont plutôt se trouver dans les peuplements éclaircis), les arbres des peuplements éclaircis ont une largeur d'aubier plus importante que ceux issus des peuplements témoins. Il ne faut pas oublier que les modèles ont été développés exclusivement à partir du site d'Ecouves. Les arbres ont tous 63 ans en 2013 et la productivité du site semble plus forte qu'ailleurs (voir section \ref{Ecouves}).  

\begin{figure}[!h]
	\centering
	\includegraphics[width=9cm]{Prediction_LA_appli_modele.pdf}	
	\caption{Prédiction des largeurs d'aubier sous le houppier pour des arbres de $D_{130}$ 30 et 50 cm (pour éviter d'utiliser le modèle en extrapolation) et pour les deux traitements "témoin" vs. "éclairci". Le graphique s'arrête à 6 m de hauteur mais la largeur d'aubier est considérée constante jusqu'à la base du houppier de chaque arbre.}
	\label{Prediction_LA_appli_modele}
\end{figure}

\FloatBarrier

\subsubsection{Modèle 4}

Le choix de développer ce modèle a été fait après le développement des modèles 2 et 3 ci-dessous et vient de l'observation de la Fig. \ref{rayon_vs_nb_cernes}. Sur cette figure il apparaissait que la quantité de duramen semblait dépendre avant tout de la taille de la rondelle. Nous avons donc cherché à développer un modèle de rayon du duramen en fonction du rayon de la rondelle de tronc considérée. La Fig. \ref{effet_satut_social_rayon_duramen_vs_rayon_rondelle} montre un effet du statut social de l'arbre et dans une moindre mesure du traitement sylvicole.

Deux formes de modèles ont été testées dont un Gompertz. Finalement nous avons retenu l'équation correspondant à la courbe affichée en bleu dans la Fig. \ref{choix_modele_4} car cette équation donnait des ajustements plus robustes aux deux extrémités de la courbe. 

La forme finale du modèle est donnée dans l'équation \ref{modele_4}.

\begin{equation} 
\boxed{	
	RayDur = a \times HBH \times (1-(b+c \times \frac{LH}{H}) \times \exp((d+e \times \exp(f \times a \times HBH)) \times RayRond))
	\label{modele_4}	
}
\end{equation} 

où $RayDur$ et $RayRond$ sont respectivement le rayon du duramen et de la rondelle en mm, $HBH$ est la hauteur de la base du houppier en m, $LH$ est la longueur du houppier en m et $H$ la hauteur totale de l'arbre en m

et avec $a = 22,31219$, $b = 0,9882738$, $c = 0,2976132$, $d = -0,001097114$, $e = -0,01549911$ et $f = -0,004962675$.


\begin{figure}[!h]
	\centering
	\includegraphics[width=9cm]{rayon_vs_nb_cernes.pdf}	
	\caption{Représentation de la quantité de duramen en fonction soit de la taille de la rondelle soit de son nombre de cernes.}
	\label{rayon_vs_nb_cernes}
\end{figure}

\begin{figure}[!h]
	\centering
	\includegraphics[width=\textwidth]{effet_satut_social_rayon_duramen_vs_rayon_rondelle.pdf}	
	\caption{Effet du statut social et du traitement sylvicole sur la relation entre le rayon du duramen et le rayon de la rondelle considérée. Dans ce graphique les rondelles pour lesquelles le rayon du duramen était nul n'ont pas été représentées.}
	\label{effet_satut_social_rayon_duramen_vs_rayon_rondelle}
\end{figure}

\begin{figure}[!h]
	\centering
	\includegraphics[width=\textwidth]{Modele4_choix_modele_suite.pdf}	
	\caption{Deux types de modèles ont été testés.}
	\label{choix_modele_4}
\end{figure}

Les prédictions du modèle pour le rayon du duramen et la largeur d'aubier sont présentées par arbre à l'Annexe \ref{rayon_dur_largeur_aubier_pred}. Les résultats apparaissent satisfaisants à ce stade. Les effets du statut social et du traitement sylvicole semblent bien pris en compte par notre modèle par le biais des variables intrinsèques à l'arbre relative au houppier qui ont été utilisées. 

Avec ce modèle, pour des arbres de même hauteur totale et à section égale de tronc, la quantité de duramen sera plus forte (et la quantité d'aubier plus forte) chez l'arbre ayant le houppier le plus court.

\FloatBarrier

\subsection{Modélisation du nombre de cernes d'aubier (modèles 2 et 3)}

Comme la largeur d'aubier est relativement constante dans le tronc sous le houppier et au-dessus de 1,30 m, et qu'il est admis que la largeur de cerne augmente du bas vers le haut de l'arbre (loi de Pressler), alors le nombre de cernes d'aubier doit logiquement diminuer entre le bas et le haut de l'arbre ou augmenter avec l'âge cambial de la section. 

\subsubsection{Modèle 2}

D'après la Fig. \ref{NbCernesDuramen_vs_NbCernesRondelle}, le nombre de cernes de duramen semble plus facile à modéliser que le nombre de cernes d'aubier. En effet, la relation entre le nombre de cernes de duramen et la hauteur dans l'arbre semble assez linéaire. La relation est encore plus linéaire si on regarde le nombre de cernes de duramen en fonction du nombre de cernes total de la rondelle. Hors cette dernière formulation peut permettre plus facilement une intégration dans un modèle de croissance du Douglas. Pour l'ajustement du modèle, nous n'avons gardé que les rondelles situées entre 1,30 m et la base du houppier.

\begin{figure}[!h]
	\centering
	\includegraphics[width=9cm]{NbCernesDuramen_aubier_vs_NbCernesRondelle.pdf}	
	\caption{Nombre de cernes de duramen et d'aubier en fonction du nombre de cernes de la rondelle à une hauteur donnée par statut social. Toutes les rondelles sont représentées, de la plus basse prise à 0,30 m à la dernière prise dans le houppier.}
	\label{NbCernesDuramen_vs_NbCernesRondelle}
\end{figure}

La relation est relativement linéaire tout le long de la tige (sauf dans le houppier à partir du moment où il n'y a plus que de l'aubier) mais nous avons fait le choix d'ajuster le modèle uniquement dans la partie du tronc située entre 1,30 m et la base du houppier (Équation \ref{NbCernesDuramen}). L'extrapolation ne posera \textit{a priori} pas de gros problème.

Les pentes sont relativement constante (Fig. \ref{NbCernesDuramen_vs_NbCernesRondelle}). Une analyse des corrélations a confirmé cette observation. La variable la plus corrélée à la pente $b$ est la largeur des 5 derniers cernes mais cette corrélation est très faible ($r=0,41$, $\text{p-value}=0,0316$). La corrélation de la pente avec le $C_{130}$ n'est pas significative.

En revanche, les variables les plus corrélées à l'ordonnée à l'origine $a$ sont : la hauteur totale de l'arbre ($r=0,63$, $\text{p-value}=0,00041$), le $C_{130}$  ($r=0,53$, $\text{p-value}=0,00448$) puis la largeur moyenne de cernes à 1,30 m ($r=0,52$, $\text{p-value}=0,00570$). Nous retiendrons le $C_{130}$ toujours pas souci de simplification.

L'effet du traitement sylvicole a ensuite été testé dans le modèle et s'est avéré être significatif sur la valeur de l'ordonnée à l'origine de la relation avec le nombre de cernes (Équation \ref{NbCernesDuramen_2}). Le $R^{2}$ de ce modèle final est de 0,94. La RMSE est de 2,16 cernes et la RMSE relative est de 9,1 \% pour les cernes de duramen et de 12,3 \% pour les cernes d'aubier (toujours entre 1,30 m et la base du houppier). En utilisant le modèle en extrapolation sous 1,30 m et dans le houppier, la RMSE est de 2,16 cernes et la RMSE relative est de 14,4 \% pour la prédiction des cernes d'aubier. 

\begin{equation} 
	NbCernesDuramen = a + b \times NbCernesTot 
	\label{NbCernesDuramen}
\end{equation}

\begin{equation} 
	NbCernesDuramen = a_{traitement} + b \times NbCernesTot + c_{traitement} \times C_{130}
	\label{NbCernesDuramen_2}
\end{equation}

On sait qu'il y a une valeur du nombre de cernes total $NbCernesTot$ en-dessous de laquelle on n'aura plus de cernes de duramen mais plus que des cernes d'aubier. C'est la valeur à laquelle le modèle du nombre de cernes de duramen (Équation \ref{NbCernesDuramen_2}) coupe l'axe des abscisses et il faut par conséquent que cette valeur soit strictement positive, c'est-à-dire que l'ordonnée à l'origine du modèle soit strictement négative, ce qui équivaut à : $a_{traitement} + c_{traitement} \times C_{130} < 0$ ou encore à $ C_{130} < -a_{traitement}/c_{traitement}$. Le modèle n'aura pas de sens si on l'applique à des arbres de plus de 60 cm de diamètre dans un peuplement témoin (ce qui doit quand même se trouver assez rarement) et à des arbres de plus de 89 cm dans des peuplements éclaircis (Tableau \ref{table_estimation_parametres_NbCernesDuramen}). Quand le nombre de cernes de duramen devient négatif dans le modèle ci-dessus (Équation \ref{NbCernesDuramen_2}), alors il faut le mettre à 0. 

Par ailleurs, le modèle permet d'estimer l'âge d'initiation du duramen : c'est la valeur du nombre de cernes de la rondelle pour laquelle le nombre de cernes de duramen vaut 0. Cet âge d'initiation est de 12 ans en moyenne et tend à diminuer quand le $C_{130}$ augmente (Fig. \ref{NbCernesDuramen_vs_NbCernesRondelle}).

\begin{equation} 
\boxed{ 
	NbCernesAub(k) = k - (a_{traitement} + b \times k + c_{traitement} \times C_{130})
	\label{NbAubier}
}
\end{equation}

avec $k=NbCernesTot$ le nombre total de cernes d'une section.

\begin{table}[h]
	\caption{Estimation des paramètres du modèle correspondant à l'équation \ref{NbCernesDuramen_2}.}
	\centering 
		\begin{tabular}{rrrrrr} 
			\hline
			Traitement & $n$ & $a$ & $b$ & $c$ & Limite max de $C_{130}$ (cm) \\
			\hline 
			Tous & 27 & -19,55381 & 0,80381 & 0,08008 & 244 \\
			\hline 
			A & 9 & -26,34712 & 0,80559 & 0,13857 & 190 \\
			B & 9 & -17,43685 & 0,80559 & 0,06628 & 263 \\ 
			C & 9 & -17,01703 & 0,80559 & 0,05712 & 297 \\
			\hline
			A & 9 & -26,34896 & 0,80564 & 0,13857 & 190 \\
			B+C & 18 & -17,18294 & 0,80564 & 0,06125 & 280 \\
			\hline
		\end{tabular}
	\label{table_estimation_parametres_NbCernesDuramen}
\end{table}

Nous avons également tenté d'ajuster un modèle à partir de variables intrinsèques à l'arbre. Un des modèles retenu (\textit{stepwise} et sélection des variables ayant un effet significatif) est celui de l'équation \ref{NbCernesDuramen_3}. Le $R^{2}$ de ce modèle est de 0,96 (0,9638). La RMSE est de 1,83 et la RMSE relative est de 9,2 \% pour les cernes de duramen et de 12,2 \% pour les cernes d'aubier (application du modèle en extrapolation sous 1,30 m et au-dessus de la base du houppier). 

Le modèle de l'équation \ref{NbCernesDuramen_4} avec un nombre réduit de variables a un $R^{2}$ de 0,96 (0,9565) également. La RMSE est de 1,98 et la RMSE relative est de 9,9 \% pour les cernes de duramen et de 13,1 \% pour les cernes d'aubier (application du modèle en extrapolation sous 1,30 m et au-dessus de la base du houppier).

\begin{multline}  
	NbCernesDuramen = a + b \times NbCernesTot + c \times Hauteur\_totale(m) + \\d \times Hauteur\_base\_houppier(m) + e \times C_{130}(cm) + \\f \times Projection\_houppier(m^{2}) + g \times LC5\_moy\_130(mm) + h \times H/D + \\i \times Longueur\_relative\_houppier 
	\label{NbCernesDuramen_3}
\end{multline} 

avec $a = -55,55465$, $b = 0,81518$, $c = -1,74565$, $d = 3,63794$, $e = -0,08807$, $f = 0,06558$, $g = 0,98904$, $h = -0,11138$ et $i = 116,3069$.

\begin{multline}  
	NbCernesDuramen = a + b \times NbCernesTot + c \times Hauteur\_totale(m) + \\d \times Hauteur\_base\_houppier(m) +  e \times LC5\_moy\_130(mm) 
	\label{NbCernesDuramen_4}
\end{multline}

avec $a = -36,62065$, $b = 0,81236$, $c = 0,37859$, $d = 0,41999$ et $e = 1,09116$. 

\FloatBarrier

\subsubsection{Modèle 3}

Ici on cherche à modéliser directement le nombre de cernes d'aubier à partir d'une formulation dynamique issue de \cite{Longuetaud2006} (Équations \ref{NBCerAubier} et \ref{NBCerAubier2}). Le nombre de cernes d'aubier à l'âge cambial $k+1$ est égal au nombre de cernes d'aubier à l'âge cambial $k$, plus le nouveau cerne qui est nécessairement un cerne d'aubier, moins une proportion des anciens cernes d'aubier qui se duraminisent. De plus si l'âge cambial est nul alors le nombre de cernes d'aubier doit être nul. Bien qu'il y ait des cernes d'aubier jusqu'à l'apex, ce qui n'est pas le cas pour le duramen, ce modèle n'est \textit{a priori} valable que dans la partie de la tige où on commence à voir de la duraminisation et où on peut considérer que le système est à l'équilibre et de plus l'hypothèse est faite que la proportion de cernes d'aubier se duraminisant est la même quels que soit l'âge cambial et l'année considérés. On ajustera donc plutôt le modèle entre 1,30 m et la base du houppier. 

\begin{equation} 
	NbCernesAub(k+1)=NbCernesAub(k)+1-NbCernesAub(k)/\alpha
	\label{NBCerAubier}
\end{equation}

\begin{equation} 
	NbCernesAub(k)=\frac{\alpha^{k}-(\alpha-1)^{k}}{\alpha^{k-1}}
	\label{NBCerAubier2}
\end{equation}

On trouve $\alpha = 20,6455$ (20,5491 si on ajuste à l'ensemble des données incluant les rondelles prises à l'empattement et dans le houppier). Ce qui signifie qu'un cerne reste pendant environ 21 ans de l'aubier avant de se duraminiser.

Après l'ajustement du modèle par arbre et analyse des corrélations, on trouve que $\alpha$ est corrélé, entre autres, à la hauteur totale ($r=-0,83$, $\text{p-value}=6,265e-08$), à la largeur moyenne des 5 derniers cernes ($r=-0,82$, $\text{p-value}=1,858e-07$), à la largeur moyenne de cerne à 1,30 m ($r=-0,75$, $\text{p-value}=5,648e-06$) et au $C_{130}$ ($r=-0,74$, $\text{p-value}=9,683e-06$). Toujours pas souci de simplifier l'utilisation du modèle, nous avons décidé de garder le $C_{130}$ comme variable explicative du paramètre $\alpha$ (Équation \ref{alpha_vs_C130} et Fig. \ref{alpha_vs_C130_figure}). Il y a de plus un effet du traitement sylvicole ("A" vs. "B+C") sur les paramètres $a$, $b$ et $c$ du modèle (utilisation de la fonction \textit{nlsList} sous R ; Fig. \ref{alpha_vs_C130_figure_effet_traitement}). La RMSE de ce modèle est de 1,87 cernes d'aubier (ou de duramen) et la RMSE relative est de 10,7 \% pour la prédiction des cernes d'aubier. Si nous extrapolons le modèle sous 1,30 m et dans le houppier, la RMSE est de 1,90 cernes et la RMSE relative est de 12,6 \%. Du point de vue des erreurs, ce modèle apparaît légèrement meilleur que le précédent.

\begin{equation} 
	 \alpha=a+\exp(b-C_{130})^c
	\label{alpha_vs_C130}
\end{equation}

\begin{figure}[!h]
	\centering
	\includegraphics[width=9cm]{alpha_vs_C130.pdf}	
	\caption{Relation entre le paramètre $\alpha$ et la variable $C_{130}$.}
	\label{alpha_vs_C130_figure}
\end{figure}

\begin{figure}[!h]
	\centering
	\includegraphics[width=9cm]{alpha_vs_C130_effet_traitement_sylvi.pdf}	
	\caption{Visualisation de l'effet du traitement sylvicole sur les paramètres $a$, $b$ et $c$ du modèle complet (Équations \ref{NBCerAubier} et \ref{alpha_vs_C130}).}
	\label{alpha_vs_C130_figure_effet_traitement}
\end{figure}

\begin{equation}
\boxed{ 
	NbCernesAub(k)=\frac{(a_{trait.}+\exp(b_{trait.}-C_{130})^c_{trait.})^{k}-(a_{trait.}+\exp(b_{trait.}-C_{130})^c_{trait.}-1)^{k}}{(a_{trait.}+\exp(b_{trait.}-C_{130})^c_{trait.})^{k-1}}
	\label{NBCerAubier3}
}
\end{equation}

\begin{table}[h]
	\caption{Estimation des paramètres du modèle correspondant à l'équation \ref{NBCerAubier3}.}
	\centering 
		\begin{tabular}{rrrrr} 
			\hline
			Traitement & $n$ & $a$ & $b$ & $c$ \\
			\hline 
			Tous & 27 & 15,82294 & 154,78830 & 0,03835 \\
			\hline 
			A & 9 & 9,70228 & 215,82210 & 0,02422 \\
			B & 9 & 16,22908 & 130,04740 & 0,05310 \\ 
			C & 9 & 18,31516 & 117,48700 & 0,09563 \\
			\hline
			B+C & 18 & 17,17231 & 131,41490 & 0,05250 \\
			\hline
		\end{tabular}
	\label{table_estimation_parametres_NbCernesAubier}
\end{table}

\FloatBarrier


\subsection{Conclusions}

\subsubsection{Évaluation des différents modèles}

La Fig. \ref{Prediction_NbCerAub_appli_modele} montre la modélisation du nombre de cernes d'aubier en fonction du nombre total de cernes à une hauteur donnée pour des arbres de diamètre à 1,30 m 30 cm et 50 cm, respectivement, issus de peuplements éclaircis ou non. 

Les Fig. \ref{Prediction_NbCerAub_appli_modele_exemple_arbres} et \ref{Prediction_LA_appli_modele_exemple_arbres} illustrent le comportement des modèles sur deux arbres contrastés (un arbre dominé du peuplement témoin et un arbre dominant du peuplement le plus éclairci) qui ont servi à ajuster le modèle (une validation en bonne et due forme sera réalisée ultérieurement sur un autre jeu de données), pour la prédiction du nombre de cernes d'aubier et de la largeur d'aubier, respectivement. Les résultats pour l'ensemble des arbres sont présentés en Annexes \ref{rayon_dur_largeur_aubier_pred}, \ref{largeur_aubier_pred} et \ref{nb_cernes_aubier_pred}, respectivement. Les modèles prédisant des nombres de cernes permettent d'obtenir un profil de largeur d'aubier (après conversion en largeur à partir des nombres de cernes et des largeurs de cernes) qui colle mieux aux données réelles avec une forme plus bombée dans le houppier et un renflement sous la base du houppier. Le modèle 1 qui prédit une largeur constante d'aubier entre 1,30 m et la base du houppier ne permet bien évidemment pas de modéliser ces variations. Le modèle 4 permet d'obtenir ces variations mais de façon parfois insuffisante.   

\begin{figure}[!h]
	\centering
	\includegraphics[width=9cm]{Prediction_NbCerAub_appli_modele_discret.pdf}	
	\caption{Prédiction du nombre de cernes d'aubier sous le houppier pour des arbres de $D_{130}$ 30 et 50 cm (pour éviter d'utiliser le modèle en extrapolation) et pour les deux traitements "témoin" vs. "éclairci" (Équations \ref{NbAubier} pour le modèle 2 et \ref{NBCerAubier3} pour le modèle 3).}
	\label{Prediction_NbCerAub_appli_modele}
\end{figure}

\begin{figure}[!h]
	\centering
	\includegraphics[width=9cm]{Prediction_NbCerAub_appli_modele_exemple_arbres.pdf}	
	\caption{Visualisation du comportement des modèles pour la prédiction du nombre de cernes d'aubier chez un arbre dominé du peuplement témoin à gauche (arbre 10) et chez un arbre dominant du peuplement le plus fortement éclairci à droite (arbre 25).}
	\label{Prediction_NbCerAub_appli_modele_exemple_arbres}
\end{figure}

\begin{figure}[!h]
	\centering
	\includegraphics[width=9cm]{Prediction_LA_appli_modele_exemple_arbres.pdf}	
	\caption{Visualisation du comportement des modèles pour la prédiction de la largeur d'aubier chez un arbre dominé du peuplement témoin à gauche (arbre 10) et chez un arbre dominant du peuplement le plus fortement éclairci à droite (arbre 25).}
	\label{Prediction_LA_appli_modele_exemple_arbres}
\end{figure}

\begin{figure}[!h]
	\centering
	\includegraphics[width=\textwidth]{comp_pred_modele2_modele3.pdf}	
	\caption{Comparaison des prédictions du nombre de cernes d'aubier pour les modèles 2 (en rouge) et 3 (en bleu). La droite tracée en noir correspond à la droite $y = x$.}
	\label{comp_pred_modele2_modele3}
\end{figure}

Les nombres de cernes d'aubier des modèles 2 et 3 sont convertis en largeurs d'aubier à partir des mesures de largeurs de cernes mesurées afin de calculer pour les trois modèles présentés une RMSE en largeur d'aubier entre 1,30 m et la base du houppier. Les RMSE des modèles 1,2, 3 et 4 sont respectivement de 5,28 mm, 5,50 mm, 5,12 mm et 5,11 mm. Les modèles 3 et 4 présentent l'erreur moyenne la plus faible et permettent une utilisation en extrapolation dans le houppier que le modèle 1 ne permet pas. 
Si on extrapole sous 1,30 m et dans le houppier on obtient pour les modèles 2, 3 et 4 des RMSE de 9,99 mm, 5,66 mm et 5,35 mm, respectivement. Le modèle 4 apparaît être notre meilleur modèle pour une prédiction des quantités d'aubier et de duramen. Les bons résultats des modèles 2 et 3, basés sur des nombres de cernes, sont aussi dûs au fait que des largeurs de cernes mesurées (et non prédites) ont été utilisées pour la conversion en largeurs.

Les différents modèles ont permis de montrer que pour des arbres de même $C130$ et de même âge, la quantité d'aubier était plus forte dans les peuplements éclaircis que dans le peuplement témoin. Nous avons également mis en évidence (modèle 4) qu'à sections de tronc égales, il y avait plus d'aubier chez un arbre dont le houppier était plus développé donc en général chez les peuplements éclaircis.

\FloatBarrier

\subsubsection{Les quantités d'aubier et de duramen - effet taille ou effet âge ?}

Quand on cherche à modéliser directement le rayon du duramen, le rayon de la rondelle explique la plus grande partie de la variabilité soit 97 \%. Le nombre de cernes de la rondelle explique quant à lui 76 \% de la variabilité quand il est introduit seul dans le modèle et n'ajoute que 2\% de $R^{2}$ quand il est utilisé en plus du rayon de la rondelle. Le nombre de cernes de la rondelle permet de différencier le haut d'un arbre qui pousse bien du bas d'un arbre qui pousse mal et qui dans notre échantillon ont le même âge mais cela n'apporte visiblement pas grand chose et il s'agit surtout d'un effet taille.


\subsubsection{Suite du travail et perspectives}

La suite de ce travail consistera à tester la validité des modèles sur d'autres jeux de données incluant notamment des arbres d'âges différents, et le cas échéant à les adapter. En effet, tous les arbres du jeu d'Ecouves ont le même âge et par conséquent le $C_{130}$ est parfaitement corrélé à la largeur de cernes moyenne. Probablement qu'il conviendrait d'inclure une variable telle que la largeur moyenne de cernes (ou des 5 derniers cernes) dans le modèle pour distinguer des arbres jeunes ayant poussé vite, d'arbres plus âgés ayant poussé plus lentement et qui pourraient avoir des $C_{130}$ du même ordre de grandeur. Cela permettrait de mieux quantifier l'effet âge qui semble pour l'instant relativement limité par rapport à l'effet taille.

Au final, quatre modèles ont été développés permettant de représenter la distribution de l'aubier et du duramen dans le tronc : un qui modélise directement la largeur d'aubier à partir du $C_{130}$ et du traitement sylvicole, un qui modélise le rayon du duramen à partir du rayon de la section et de variables hauteur de la base du houppier et longueur relative du houppier et deux qui prennent en entrée le nombre total de cernes de la section en plus du $C_{130}$ et du traitement sylvicole. Il est possible qu'un de ces modèles se comporte mieux qu'un autre en extrapolation sur un nouveau jeu de données.

Nous avons fait le choix d'utiliser le $C_{130}$ dans la plupart de nos modèles pour qu'ils soient plus facilement applicables. Cependant l'utilisation d'autres variables explicatives comme par exemple la largeur moyenne à 1,30 m des 5 derniers cernes pourrait permettre d'améliorer encore les résultats. Quelques autres modèles sont proposés en ce sens à titre d'illustration.

Enfin nous avons choisi de considérer dans une partie des modèles deux traitements sylvicoles : "témoin" vs. "éclairci" ce qui n'est pas satisfaisant. Nous devons réfléchir à une variable permettant de mieux quantifier l'intensité de l'éclaircie, peut-être basée sur des trajectoires de $RDI$ ou sur les distances aux voisins. Par ailleurs, des modèles sont également proposés qui ne prennent pas en entrée la variable "traitement sylvicole" mais uniquement des variables intrinsèques à l'arbre. 

Les modèles pourront également être raffinés avec l'utilisation de modèles mixtes pour la prise en compte de l'effet aléatoire "arbre". 

Une fois les modèles jugés suffisamment robustes, ils pourront être intégrés dans le simulateur de croissance \textit{SimCop} sous \textit{Capsis}.

Une poursuite nécessaire du travail serait de rechercher des jeux de données à l'étranger (notamment Allemagne et Belgique) incluant des mesures de largeur d'aubier à plusieurs hauteurs mais également des mesures à 1,30 m mais sur un grand nombre d'arbres représentant des conditions de croissance variées.


%\bibliographystyle{elsarticle-harv}
%\bibliographystyle{frplain}
\bibliographystyle{francaissc}
\bibliography{bibliographie}

\newpage

\appendix

\chapter{Revue des équations de profils de tiges extrait de \cite{Li2010}} 

\begin{sidewaysfigure}[!h]
	\begin{center}
	\includegraphics[width=1.\textwidth]{Li_Weiskittel_2010_Table2}
	\end{center}
	\caption{Extrait de \cite{Li2010}}
	\label{Li_Weiskittel_2010_Table2}
\end{sidewaysfigure}

\chapter{Variations de $\theta$ en fonction de la hauteur relative} \label{AppendixA}
\includepdf[pages = {1-3}]{theta_vs_hauteur_relat.pdf}

\chapter{Accroissement ligneux, ajustement de la fonction $ShX$ et calcul des dérivées 1\iere, 2\ieme et 3\ieme} \label{Appendix_Shx}

Les indicateurs $H_f$, $H_s$ et $H_b$ sont définis de la façon suivante :
\begin{itemize}
\item $H_s$ est tout d'abord calculé à partir de la dérivée troisième : c'est la hauteur correspondant au premier point en partant de l'apex pour lequel la dérivée troisième s'annule, ou à défaut au point pour lequel elle est maximale ;
\item  $H_f$ est ensuite calculé à partir de la dérivée seconde : c'est la hauteur correspondant au premier point en partant de l'apex pour lequel la dérivée seconde s'annule, ou à défaut à la valeur de $H_s$ ;
\item $H_b$ est enfin calculé à partir de la dérivée première : par construction, la dérivée est nulle au niveau de l'apex ; $H_b$ est la hauteur correspondant au premier point suivant pour lequel la dérivée s'annule. A défaut, $H_b$ est fixé à la hauteur correspondant au second point pour lequel la dérivée seconde s'annule, ou à 0 si la dérivée seconde ne s'annule qu'une fois, ou enfin à la valeur de $H_f = H_s$ si elle ne s'annule jamais.
\end{itemize}

\section{Site d'Ecouves pour l'année 2010}
\includepdf[nup=2x3, pages = {1-27}]{ShX_et_derivees_par_arbre_Ecouves_2010.pdf}
\section{Site d'Amance pour l'année 1991}
\includepdf[nup=2x3, pages = {1-24}]{ShX_et_derivees_par_arbre_Amance_1991.pdf}

\chapter{Accroissement ligneux et évolution de l'ajustement de la fonction $ShX$ au cours du temps} \label{Appendix_Shx_2}
\section{Site d'Ecouves : ajustements pour les années 1960, 1970, 1980, 1990, 2000 et 2010}
\includepdf[pages = {1-3}]{profils_acc_par_arbre_comparaison_modele_Ecouves.pdf}
\section{Site d'Amance : ajustements pour les années 1965, 1970, 1975, 1980, 1985 et 1990}
\includepdf[pages = {1-4}]{profils_acc_par_arbre_comparaison_modele_Amance.pdf}

\chapter{Remontées des différentes hauteurs : hauteur totale, Hf, Hs, Hb et Hf2} 
La hauteur de la base du houppier prédite rétrospectivement à partir de l'ajustement des profils est représentée en orange.
\label{Appendix_Remontees}
\section{Site d'Ecouves}
\includepdf[pages = {1-3}]{Remontees_H_Hf_Hs_Hb_par_traitement_par_statut.pdf}
\section{Site d'Amance}
\includepdf[pages = {1-4}]{Remontees_H_Hf_Hs_Hb_par_traitement_par_statut_Amance.pdf}

\chapter{Résidus du modèle d'incrément en surface de cerne par site, par arbre et par année calendaire} 
\label{Residus_modele_increment_par_arbre_par_annee}
\section{Site d'Ecouves}
\includepdf[pages = {1-3}]{residus_par_annee_par_arbre_Ecouves.pdf}
\section{Site d'Amance}
\includepdf[pages = {1-4}]{residus_par_annee_par_arbre_Amance.pdf}

\chapter{Ajustement du modèle prédictif aux données} 
\label{Ajustement_modele_predictif}
\section{Site d'Ecouves : ajustements pour les années 1960, 1970, 1980, 1990, 2000 et 2010}
\includepdf[pages = {1-3}]{profils_acc_par_arbre_comparaison_modele_Ecouves_predictif_modele_7_param_mixte.pdf}
\section{Site d'Amance : ajustements pour les années 1965, 1970, 1975, 1980, 1985 et 1990}
\includepdf[pages = {1-4}]{profils_acc_par_arbre_comparaison_modele_Amance_predictif_modele_7_param_mixte.pdf}

\chapter{Evolution avec l'âge des profils normés de surface de cerne prédits par la partie fixe du modèle pour les arbres mesurés.} 
\label{Annexe_profils_normes}
\section{Site d'Ecouves}
\includepdf[pages = {1-3}]{profils_acc_par_arbre_comparaison_modele_mixte_arbre_partie_fixe_2_ecouves}
\section{Site d'Amance}
\includepdf[pages = {1-4}]{profils_acc_par_arbre_comparaison_modele_mixte_arbre_partie_fixe_2_amance}

\chapter{Rayon du duramen et largeur d'aubier} 
\label{Appendix_LA}
\includepdf[pages = {1-3}]{profils_aubier_largeur.pdf}

Les lignes grises en pointillés représentent dans l'ordre : la première branche verte, le premier verticille vert (3/4 des branches vertes) et la hauteur totale de l'arbre. La ligne verte indique l'endroit où la largeur de l'aubier devient supérieure au rayon du duramen.

\chapter{Surfaces du duramen et de l'aubier} 
\label{Appendix_surface_aubier}
\includepdf[pages = {1-3}]{profils_aubier_surface.pdf}

Les lignes grises en pointillés représentent dans l'ordre : la première branche verte, le premier verticille vert (3/4 des branches vertes) et la hauteur totale de l'arbre. La ligne verte indique l'endroit où la surface de l'aubier devient supérieure à celle du duramen.

\chapter{Nombre de cernes du duramen et de l'aubier} 
\label{Appendix_nb_cernes_aubier}
\includepdf[pages = {1-3}]{profils_aubier_cernes.pdf}

Les lignes grises en pointillés représentent dans l'ordre : la première branche verte, le premier verticille vert (3/4 des branches vertes) et la hauteur totale de l'arbre. La ligne verte indique l'endroit où le nombre de cernes de l'aubier devient supérieure à celui du duramen.

\chapter{Proportions en largeur, en surface et en nombre de cernes de l'aubier} 
\label{Appendix_prop_aubier}
\includepdf[pages = {1-3}]{profils_aubier_proportions.pdf}

Les lignes grises en pointillés représentent dans l'ordre : la première branche verte, le premier verticille vert (3/4 des branches vertes) et la hauteur totale de l'arbre. La ligne verte à 0.5 indique l'endroit où la proportion d'aubier devient supérieure à celle de duramen pour chacune des trois variables étudiées.

\chapter{Évolution de la largeur d'aubier à l'empattement} 
\label{largeur_aubier_empattement}
\includepdf[pages = {1-3}]{profils_aubier_largeur_empattement.pdf}

\chapter{Analyse des résidus du modèle permettant de prédire la largeur moyenne d'aubier entre 1,30 m et la base du houppier} 
\label{largeur_aubier_tronc_residus}

\begin{figure}[!h]
	\centering
	\includegraphics[width=9cm]{LA_vs_C130_tronc_residus_vs_traitement.pdf}	
	\caption{Résidus du modèle (Équation \ref{modele_2_LA_sous_houppier_sans_traitement}) en fonction du facteur traitement sylvicole (A = peuplement témoin ; B = traitement intermédiaire ; C = éclaircies intensives).}
	\label{LA_moy_residus_vs_traitement}
\end{figure}

\begin{figure}[!h]
	\centering
	\includegraphics[width=\textwidth]{LA_vs_C130_tronc_residus_vs_variables_arbre.pdf}	
	\caption{Résidus du modèle (Équation \ref{modele_2_LA_sous_houppier_sans_traitement}) en fonction des variables au niveau "arbre".}
	\label{LA_moy_residus_vs_variables_arbre}
\end{figure}

\chapter{Évolution de la largeur d'aubier dans le houppier} 
\label{largeur_aubier_houppier}
\includepdf[pages = {1-3}]{profils_aubier_largeur_houppier.pdf}

\chapter{Prédictions du rayon du duramen et de la largeur d'aubier à partir du modèle 4} 
\label{rayon_dur_largeur_aubier_pred}
\includepdf[pages = {1-3}]{profils_RayDur_par_arbre_modele_final_2.pdf}

\chapter{Visualisation du comportement des modèles 1, 2, 3 et 4 pour la prédiction de la largeur d'aubier} 
\label{largeur_aubier_pred}
\includepdf[pages = {1-3}]{largeur_aubier_appli_modeles_tous_arbres.pdf}

\chapter{Visualisation du comportement des modèles 2 et 3 pour la prédiction du nombre de cernes d'aubier} 
\label{nb_cernes_aubier_pred}
\includepdf[pages = {1-3}]{nb_cernes_aubier_vs_nbcernes_appli_modeles_tous_arbres.pdf}


\end{document}


